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ABSTRACT 

The United States Coast Guard (USCG) is in the process of upgrading the 
hardware of the Loran-C Radionavigation System Control System. As part of this 
effort, the Computer- Assisted Loran-C Controller (CALOC), is also in need of im- 
provement. CALOC performs four tasks: abnormality detection, time difference 
control, recordkeeping, and blink control. The work reported in this thesis focuses on 
time difference control. In many instances, CALOC does not accurately control the 
time difference error (TDE) within the established USCG control procedures. Two 
new algorithms are proposed here to control TDE more effectively: a proportional- 
integral-derivative (PID) controller and a Kalman filter. Actual TDE data recorded 
at three different master stations covering five Loran-C chains is used to evaluate the 
performance of the proposed controllers. The PID controller shows a substantial im- 
provement in control compared to CALOC, and the Kalman filter exhibits even better 
performance, based on preliminary results. This improvement in control correlates 
directly with an increase in both predictable accuracy and repeatable accuracy. 
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I. 



INTRODUCTION 



As part of the 1996 Federal Radionavigation Plan, Loran-C is scheduled to be 
decommissioned on 31 December 2000, and the Global Positioning System (GPS) is 
to become the sole source of global radionavigation. However, there is still significant 
interest in the use of Loran-C for air, land and marine navigation, and additional 
concern about not having a back-up, or complementary, navigation system to GPS. 
Under the Coast Guard Authorization Act of 1996, the Department of Transportation 
is required to present a report on the impact of the Loran-C decommissioning. There 
are numerous factors that will influence the decommissioning decision, including in- 
ternational treaties and alliances, costs of maintaining more than one federally-funded 
radionavigation system, and domestic user requirements. Improving the accuracy of 
Loran-C and reducing the manpower required to operate the stations will lead to re- 
duced operational costs, thus making Loran-C a feasible option to complement GPS. 

Maintaining the technological edge of legacy systems has always been a chal- 
lenge to design engineers. The United States Coast Guard’s (USCG) efforts to im- 
prove the timing control of the Loran-C transmitted pulse is no different. In the 
late 1980s, the USCG embarked upon an aggressive multi-year project, titled “Elec- 
tronics Equipment Replacement Project” (EERP), to upgrade the entire Loran-C 
system. As part of EERP, Coast Guard engineers at the Electronics Engineering 
Center (EECEN) have identified commercial-off-the-shelf (COTS) items that will re- 
place outdated control subsystem hardware and consolidate personnel tasking. 

To complement this hardware upgrade, the control algorithm needs to be im- 
proved. The existing system, the Calculator- Assisted Loran-C Controller (CALOC), 
does not maintain Loran-C within USCG control specifications. For example, in the 
recorded data provided by the USCG for use in this thesis, there are cumulative time 
difference errors (TDEs) in excess of 100 ns and instantaneous TDE exceeding the 
USCG established control limit of ±50 ns. These large TDEs are directly related 
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to the inadequacy of the existing control algorithm. The Loran-C operators, know- 
ing that CALOC is unreliable, do not allow it to automatically control the Loran-C 
chain. It is impossible to determine how the operators’ actions influence the accuracy 
of CALOC. 

Timing control of the Loran-C transmitted pulse is directly related to both 
the absolute accuracy and the repeatable accuracy of the Loran-C system. Both 
types of accuracy are important, but for different reasons. For example, the absolute 
accuracy is important to the mariner who may be entering an unfamiliar port for 
the first time while the repeatable accuracy is important to the coastal fisherman 
who lays lobster traps in the morning and wants to return to the same location in the 
evening. In general terms, the absolute accuracy of the Loran-C system includes both 
the random errors and the systematic errors while the repeatable accuracy includes 
only the random errors. The algorithm for control of the transmitted pulse timing 
is only capable of controlling the random errors. Hence, an improvement in control 
of the transmitted pulse timing leads directly to an improvement in the repeatable 
accuracy. 

The goal of this work is to develop a new control algorithm that will maintain 
control of both the instantaneous and cumulative TDE within USCG specifications, 
and in turn improve operator confidence in the control system. The algorithm is 
developed using actual data received from the monitor stations and the output from 
the existing control algorithm. The data is first preprocessed to remove as much of 
CALOC’s influence as possible to achieve uncorrected data. The preprocessed data 
retains the effects of atmospheric, terrain-related and man-made noise that cause 
random errors in the Loran-C signal. This simplifies the model required for imple- 
menting the algorithm, and it also enables a more accurate solution since the control 
algorithm is tested using field-recorded data. 

In order to develop a new control algorithm, an understanding of Loran-C 
with an emphasis on CALOC is required. A brief overview of Loran-C is contained in 
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Chapter II, along with a description of the EERP as it applies to Control Equipment 
Consolidation. Chapter III contains pertinent information about CALOC, detailing 
its algorithm and the theory supporting it. In Chapter IV, a proportional-integral- 
derivative (PID) control algorithm is developed which leads to a data-dependent 
model discussed in Chapter V. Finally, conclusions and recommendations for further 
study are given in Chapter VI. 

The appendices contain Matlab plots and the code used to generate them. The 
strip charts based on the recorded data from CALOC are contained in Appendix A. 
Appendix B contains the plots for the same data controlled by the PID controller, 
and Appendix C contains the plots for the Kalman filter. Finally, the Matlab code is 
included in Appendix D, beginning with a brief description of the code and definitions 
of functions and variables. 
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II. 



OVERVIEW OF LORAN-C 



LORAN, an acronym for LOng RAnge Navigation, is a member of the family 
of radionavigation systems that provides geographic position-fixing data in the form 
of hyperbolic lines of position. It operates on the principle that radio waves travel 
at a constant velocity. The measurement of the difference in arrival time of radio 
pulses from two time-synchronized transmitting stations establishes a hyperbolic line 
of position. The intersection of two or more lines of position provides a fix, or location. 
[Ref. 1] 

A. HISTORY 

Loran-C can trace its roots back to a hyperbolic radionavigation system pro- 
posed by R. J. Dippy in 1937 and later implemented as the British Gee system in 
early 1942. During this period in the United States, the same principles of hyperbolic 
radionavigation were being researched by the MIT Radiation Laboratory. Spurred by 
the events of World War II, development of Loran proceeded rapidly in the United 
States with a chain of transmitters, later called standard Loran or Loran- A, that were 
operated by the USCG starting in January 1943. These initial stations, along with 
others operated by the Royal Navy and the Royal Canadian Navy, provided coverage 
over the North Atlantic convoy routes. Additional stations were constructed in the 
Aleutian Islands and central Pacific Ocean so that by the end of the war Allied mili- 
tary ships and aircraft enjoyed nighttime coverage of over 30% of the earth’s surface. 
[Ref. 1, 5] 

Following World War II, research sponsored by the Department of Defense 
(DoD) and other agencies explored a more accurate and longer range version of Lo- 
ran. This research led to various improvements, including low frequency (LF) Loran, 
Loran-B, CYCLAN, and CYTAC. The Cycle Tactical (CYTAC) Bombing and Nav- 
igation System, consisting of a master and two slave transmitter stations, was the 
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most promising development to evolve. The CYTAC system allowed a single mobile 
receiver to determine its position by automatically measuring the difference in time 
of arrival of radio- frequency (RF) pulses from the three transmitters. [Ref. 5] 

In 1957, the United States Navy identified a need for a long-range maritime 
radio navigation system for use with the ballistic missile submarines. However, the 
accuracy and range required by the Navy considerably exceeded the capabilities of 
the existing Loran-A equipment. Based on the results of the CYTAC testing, it was 
believed that the requirement for long-range navigation could be met by implement- 
ing the CYTAC concept. CYTAC was made operational in 1957 and placed under 
the control of the USCG, designating it Loran-C, in 1958. One of the significant 
improvements was using a lower frequency band of 90 to 110 kHz, as opposed to 
the 1.85 to 1.95 MHz of Loran-A, which provided a significant increase in range for 
Loran-C. Additional technical improvements also enabled more accurate geographic 
positioning. [Ref. 1, 5] 

Today, Loran-C is a series of transmitter chains throughout the world, some 
operated by the USCG and others by host nations that provide highly accurate, 
24-hour-a-day, all-weather radionavigation used for ocean and coastal navigation. In 
the United States it is the federally provided radionavigation system for the U.S. 
Coastal Confluence Zone (CCZ).^ Loran-C is used primarily by mariners, both do- 
mestic and international, but is also widely used by civilian aircraft in the United 
States. [Ref. 1] 

B. LORAN-C SYSTEM 

The basic element of the Loran-C system is the loran chain. The chain consists 
of one master loran station and at least two secondary stations. Typical configura- 
tions for Loran-C chains are illustrated in Figure 1. The secondary stations are 

'The CCZ is defined as the area seaward of a harbor entrance to 50 nautical miles offshore or 
the edge of the Continental Shelf (100 fathom curve), whichever is greater. [Ref. 1] 
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Figure 1. Typical Configurations for Loran-C Chains. [Ref. 1] 

generally labelled starting with Zulu and working backwards through the alphabet. 
Each Loran-C chain provides signals suitable for accurate navigation over a coverage 
area. These coverage areas overlap in many parts of the United States and its coastal 
waters, where two or more chains may be received and used for navigation. [Ref. 1] 
The concept of Loran-C is based on the time difference (TD) between the 
time of arrival (TOA) of the transmitted pulses of a master station and one of several 
secondary stations in a loran chain. Beginning with the master station, each loran 
station (LORSTA) within a chain transmits a series of phase-coded pulses. The timing 
of pulse transmission allows for each station’s signal to propagate throughout the 
area of coverage before the next station transmits. Normally, the secondary stations 
transmit in alphabetical order (i.e., Victor, Whiskey, Xray, Yankee, and then Zulu). 
A receiver within the chain’s area of coverage measures and displays the elapsed time 
between the signals from each station. The TD between the master and secondary 
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stations yields a hyperbolic line of position on which the receiver must lie. A fix^ is 
determined by two or more intersecting hyperbolic lines of position. [Ref. 6] 

The chain’s master and secondary stations continue to sequentially transmit 
the phase-coded pulses with the master station transmitting a fixed time following the 
final secondary station in the chain. The length of time between each master station 
transmission is determined by the group repetition interval (GRI), which is assigned 
to each loran chain upon installation. Assigned GRIs range from 59,300 to 99,990 
microseconds, which means that a chain’s complete transmission cycle will repeat 
anywhere from 10 to 17 times each second. In order to ensure proper transmission 
timing, progressively longer emission delays^ are assigned, so the transmission of the 
phase-coded pulses from each secondary station are allowed to propagate throughout 
the chain’s area of coverage. [Ref. 6] 

A typical transmission pulse pattern for a Loran-C chain, made up of a master 
station and three secondary stations, is shown in Figure 2. The emission delay of each 
secondary station, depicted as TDX, TDY and TDZ, is the length of time between 
the start of transmission for the master station and that of each secondary station. 
The GRI is illustrated as the time between the start of transmission for the master 
station and the next time it transmits. Also depicted in Figure 2 is an exploded 
view of the Loran-C pulse shape, consisting of sine waves enclosed within a teardrop 
shaped envelope, referred to as a t-squared pulse. [Ref. 1, 6] 

C. THE ELECTRONIC EQUIPMENT REPLACEMENT 
PROJECT 

In 1989, the United States Coast Guard (USCG) embarked on an effort to 
redesign portions of the Loran-C system aimed at taking it into the 21®*' century. The 

fix is a known position determined the intersection of two or more lines of position (LOPs) 
determined from terrestrial, electronic and/or celestial data. [Ref. 1] 

^Emission delay (ED) is the time difference, in /x seconds, between a master loran station trans- 
mission and a given secondary station transmission. [Ref. 1] 
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Figure 2. Pulse Pattern and Detailed Pulse Shape for Loran-C Transmission. [Ref. 1] 



9 



USCG’s Electronics Engineering Center (EECEN) was tasked to perform a complete 
analysis of several subsystems under a multi-year project known as the Electronic 
Equipment Replacement Project (EERP) using the following criteria: the supporta- 
bility of present day equipment; the support and maintenance of existing equipment 
projected ahead five years; the desire to enhance and expand automation; the need 
to respond to new system requirements; and the desire to remain in close step with 
technology. The primary purpose of EERP is to identify immediate support prob- 
lems, develop near term solutions to these problems and chart the course for system 
improvements while addressing new requirements. Redesign of various portions of 
the Loran-C system is necessary to meet the demands of maintaining and operating 
the system into the next century. This thesis is concerned with analyzing the existing 
Loran-C Control Subsystem and developing an updated control algorithm to com- 
plement the hardware improvements being installed in support of the EERP’s Loran 
Consolidated Control Station (LCCS). [Ref. 7] 

1. Plan V - Control Equipment Consolidation 

One of the goals of the EERP is to address a major flaw in the current control 
subsystem: excessive dependence on human watchstanders to evaluate incoming data 
from multiple sources. Currently, watchstanders constantly receive and must then 
understand both visual and audio outputs from four distinct pieces of equipment. 
Watchstanders have no method to interpret the received information unless they can 
correlate the data with information from another piece of equipment. Figure 3 il- 
lustrates the numerous inputs required for the watchstanders to interpret. While 
reducing the need to manually correlate information is not an adequate solution, all 
the information should be consolidated into one piece of equipment. This consol- 
idation would significantly simplify the watchstander’s job. Figure 4 outlines the 
various modules of the control subsystem. As part of the redesign, the input/output 
relationships shown in the figure would remain unchanged. [Ref. 7] 
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2. LORAN-C Consolidated Control Station (LCCS) 

The requirements analysis associated with the Control Equipment Consolida- 
tion concluded that the watchstander is presented with too many sources of Loran 
chain status information and that this information must be consolidated into a single 
view of the Loran chain. The purpose is to generate this single view of the Loran 
environment and display it for the watchstander. The outgrowth of this project is 
the Loran Consolidated Control Station (LCCS). 

The primary mission of the LCCS is to monitor and control LORSTAs and 
Primary Chain Monitor Sets (PCMS) currently functioning as part of the Loran-C 
Radionavigation System. The secondary mission is to replace the current suite of 
chain control equipment, which includes the Remote Site Operating Set (RSOS), 
Calculator-Assisted Loran-C Controller (CALOC), Chain Recorder Set (CRS) and 
Teletype (TTY). The chain control equipment at the three Continental United States 
(CONUS) Loran Control Stations in Malone, Florida, Seneca, New York, and Mid- 
dletown, California, will be replaced. The LCCS will reduce and automate the data 
collection and decision-making functions at these stations, so they may be remotely 
controlled from central control sites in Alexandria, Virginia, and Petaluma, Califor- 
nia. [Ref. 2] 

The hardware and software components of the LCCS are divided into sev- 
eral Hardware Configuration Items (HWCI) and Computer Software Configuration 
Items (CSCI). Figure 5 depicts a top-level LCCS block diagram illustrating these 
configuration items. The LCCS Tactical Advanced Computer (TAC-4) is an HP 9000 
workstation with an HP-UX (Unix-based) operating system. Two CSCIs will be used: 
one is the executable software to run LCCS, and the other is the network manage- 
ment software. Other HWCIs are the Data Retrieval/Monitoring System and the 
Commnnications Interface. [Ref. 2] 

In concert with the implementation of Control Equipment Consolidation and 
installation of the LCCS, the LORSTAs will make a significant change in the method 
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Figure 5. LCCS Top-Level Block Diagram. [Ref. 2] 

of communicating betAveen stations. The current process involves using either dedi- 
cated or leased telephone lines. With the upgrade, communications will employ X.25 
data lines and packet assembling/disassembling equipment at all CONUS LORSTAs 
and PCMS sites. In addition to providing potentially significant cost savings, the 
packet switching system is able to be configured from any location. This will allow 
for monitoring and troubleshooting of the Loran-C communications network from a 
single location, such as EECEN. [Ref. 2] 

3. Proposed Effort 

One of the areas of the Control Subsystem that is not currently scheduled to 
be upgraded is the control algorithm resident in CALOC. The Coast Guard intends 
to port the current CALOC algorithm to the new Unix environment of LCCS with 
minimal functional change. The primary goal of this thesis is to investigate alter- 
nate solutions that will replace the existing algorithm. The new control algorithm 
must keep both the instantaneous time difference error (TDE) within ±50 ns and 
the cumulative TDE close to the controlling standard time difference (CSTD) while 
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providing the minimum number of local phase adjustments (LPA), which are trans- 
mitter timing adjustments issued in multiples of 20 ns. The basic control law that 
will be followed is to minimize the TDE by maintaining the time difference as close as 
possible to its control number while minimizing the magnitude and number of timing 
adjustments applied. [Ref. 8] 

In this chapter, a brief history and system description of Loran-C was pre- 
sented. It also outlined the EERP, focusing primarily on Control Equipment Con- 
solidation, and in particular the LCCS. As part of EERP, a new control algorithm 
needs to be developed to replace CALOC. In order to develop a suitable algorithm, 
an understanding of the existing system is required. The next chapter discusses the 
functionality of CALOC, including the four general tasks it performs. It also describes 
the control policies and the role of monitor stations in the performance of Loran-C. 
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III. CALCULATOR ASSISTED LORAN-C 

CONTROLLER 

The Calculator-Assisted Loran-C Controller (CALOC) consists of a control 
algorithm implemented as a desktop calculator and various hardware peripherals, 
which may be installed at a Loran-C time difference control station. Its function is to 
relieve the operator of many tedious and routine jobs and perform several computa- 
tions. CALOC monitors and controls the phase timing of the Loran-C signal on up to 
four baselines while the envelope timing is presumed to be monitored and maintained 
at each individual transmitting station via the installed local signal envelope-to-cycle 
difference (ECD) measurement, limit check, and alarm circuits. [Ref. 8] 

A. TASKS 

In order to monitor and control the Loran-C system, CALOC has been de- 
signed to perform four general tasks: abnormality detection; time difference (TD) 
control; recordkeeping; and blink control. In performing these tasks, CALOC acts 
mainly to assist, not replace, the human operator, especially in the area of subjective 
judgment which must be applied during many system casualty conditions. [Ref. 8] 

1. Abnormality Detection 

An equipment casualty at either the master or the secondary transmitting 
station can cause a sudden excursion, or abnormal behavior, in the received phase 
timer difference (TD). Since CALOC cannot possibly anticipate the symptoms of, 
and corrections for, all possible system casualties, it simply inhibits control task 
operation on a baseline when it finds that the baseline is abnormal. At the same 
time, CALOC notifies the human operator of the abnormal condition and then relies 
upon the operator to correct it. When the operator is satisfied that the baseline is 
again normal, he returns control to CALOC. [Ref. 5, 8] 



15 



To determine if a baseline is abnormal, CALOC uses a statistical algorithm. 
The CALOC system, in the face of continuously changing received noise power, con- 
tinually evaluates the two hypotheses — baseline normal or baseline abnormal — 
while attempting to satisfy the conflicting goals of minimizing the time to detect 
and the probability of making an error. In this manner, CALOC relieves the human 
operator of the need for constant attention to his receiver displays and strip charts. 
[Ref. 5, 8] 

2. Time Difference (TD) Control 

Under steady-state or normal conditions, the time difference will fluctuate 
due to the transmitter oscillator activity and other transmitter station variations, 
and propagation variations. Because the received signal is contaminated by noise, 
the optimum number and size of LPAs necessary to keep the time difference at the 
established control number is not obvious, and it is difficult for a human operator to 
properly weigh all of the many factors in the decision to insert an LPA. [Ref. 5, 8] 

The most important function of CALOC is to quantify the process of observing 
the present and past time difference behavior, weigh the competing factors involved 
in the insertion of an LPA, and choose the optimum correction to insert. This system 
transmits the optimum LPA value to the appropriate Remote Control Interface (RCI) 
for insertion at the remote station when operating in the automatic control mode. 
In the semi-automatic control mode, the system recommends that the operator make 
the necessary timing correction. The operator is allowed the option of inserting the 
LPA manually, allowing CALOC to insert the LPA, or aborting the recommendation 
altogether. In this way, the operator has the final decision on the LPA when CALOC 
is placed in the semi-automatic control mode. If the CALOC system is initialized 
for a default manual mode, recommendations will be made, but CALOC will not 
transmit the RCI control messages. The manual mode is provided for those control 
stations that have unreliable teletype communications. [Ref. 5, 8] 



16 



3. Recordkeeping 

In the execution of the above tasks, CALOC generates many numerical values 
which are needed by the Coordinator of Chain Operations (COCO) for the proper 
management of the Loran-C system. Since these values can assist the TD control 
station personnel in the preparation of their operational messages, CALOC writes 
them onto a printer tape and drives a plotter to show exactly how the baselines 
under its control performed. By using the printer tape and plotter output, which 
CALOC annotates with the day’s major events, the operator is relieved of much of 
the recordkeeping burden. [Ref. 5, 8] 

4. Blink Control 

Blink codes are used to transmit the status of Loran-C system accuracy and 
the reliability of the stations within the chain. If one of the stations is out of tolerance 
or unusable for some reason, the affected station will blink. Blinking is a repetitive 
on-off pattern of the first two pulses of the secondary signal, which indicates that the 
baseline is unusable for one of the following reasons: 

• Time Difference out of tolerance, 

• Envelope-to-Cycle Difference out of tolerance, 

• Improper phase code or Group Repetition Interval, or 

• Master or secondary station operating at less than one-half of specified output 

power, or master station is off-air. 

When any of the above conditions occur, the operator will initiate Blink via CALOC. 
[Ref. 5] 

B. CONTROL POLICIES 

The Coast Guard’s stated objectives in controlling the TD of a Loran-C base- 
line at the primary control monitor are to use a limited number of LPAs while keeping 
the TD close to its assigned number as well as keeping the cumulative TDE close to 
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zero. Every 15 minutes, CALOC determines whether an LPA is required to correct 
the time difference. To convert the evaluation or determination of whether or not 
to insert an LPA into a numerical procedure suitable for CALOC, a control cost is 
defined. This cost, J, is a weighted sum of the deviation of the time difference from 
the assigned number, the deviation of the cumulative TD error from zero, and the 
size of the LPA being considered: 

{ZN-Nc + uf e^u) 

Kf (20ns)2 

where Zj^r is the most recent TD estimate, u is the value of the individual LPA being 
evaluated, e{u) is the total predicted cumulative error, {Zj\r — Nc) is our predicted 
TDE, and K{ and Kg are parameters used to alter the controller function. The cost 
is computed using a least squares method for each possible LPA. The LPA which 
results in the least cost is then selected. [Ref. 5, 8] 

In order to assist in minimizing the number of LPAs during times of high at- 
mospheric noise, Loran uses two different control policies (CONPOL). Two sets of 
weighting factors, Ki and Kg, associated with CONPOL 1 and CONPOL 2, respec- 
tively, are stored in CALOC and used to provide adaptive properties to the algorithm. 
The time of their application is controlled by the COCO. One use of the control poli- 
cies is to provide close control during the day and loose control at night when there 
is, in general, more noise in the atmosphere. [Ref. 5, 8] 

C. MONITORING STATIONS 

The purpose of an unmanned remote Loran-C monitor site is to receive signals 
from Loran-C transmitting stations and provide regular reports on cycle time dif- 
ference, absolute envelop- to-cycle difference (ECD), signal-to-interference ratio, and 
signal strength to the respective chain control sites. The chain control sites use the 
monitor information to hold Loran-C transmitting station signals to within assigned 
tolerances. Each Loran-C monitor site consists of a receiving whip antenna with a 
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wire ground plane and a small, all-weather shelter, which houses the Loran-C monitor 
equipment and attendant prime power and communication service tie points. There 
are two primary pieces of equipment that the monitoring station uses to process the 
received Loran-C signals: the Austron 5000 and the PDP-8/E computer. 

1. Austron 5000 

The Austron 5000 is a high accuracy receiver used as a remote sensor for 
control of the Loran-C chain. The receiver may be used on a time-shared basis to 
monitor and track up to four chains consisting of eight stations. The system consists 
of a receiver section which filters and amplifies the Loran-C signal and then digitizes 
samples upon command from the computer section. The receiver section is initialized 
and controlled by a PDP-8/E general purpose 12-bit minicomputer. All processing is 
done in digital servo loops and filters written in the monitor program. [Ref. 5] 

The receiver uses direct RE sampling. During track mode, three samples are 
taken on each pulse: one zero crossing sample at 30 /rs called the phase strobe; one 
at -1-2.5 fis from the phase strobe called the amplitude strobe; and one at -7.5 //s 
from the phase strobe called the cycle strobe. Additionally, a guard sample is taken 
from the pulse on a time-shared basis to check for early energy. These three sample 
and hold circuits are followed by a 12-bit analog-to-digital (A/D) converter which 
sequentially digitizes the three samples and passes them to the PDP-8/E computer 
for processing in software phase-locked loops and filters. [Ref. 5] 

2. PDP-8/E Computer 

Each monitor station has two PDP-8/E computers: one is used to control the 
Austron 5000, and the other is maintained in a standby condition. The in-service 
computer treats the receiver section, including the Austron 5000, as a peripheral 
device and collects TD data, envelope data, receiver gain and noise data from the 
Austron 5000. This information is formatted by the computer into the data messages. 
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coded strip chart update messages, and error messages. These messages are then sent 
to the control monitor station via a patch panel and modem. [Ref. 5] 

D. TIME DIFFERENCE ERROR (TDE) 

The time difference error (TDE) is the principal parameter used to control the 
timing of the transmitted Loran-C pulse. There are two components that make up the 
TDE; the controlling standard time difference (CSTD) and the time difference (TD). 
The CSTD is established through a chain calibration process, of which the primary 
purpose is to ensure that the emission delay (ED) of each secondary station is set to 
the value published by the USCG. The time difference is a measurement taken at the 
monitor station and is used to determine the time difference error (TDE). [Ref. 9] 

1. Controlling Standard Time Difference (CSTD) 

The CSTD for each master-secondary pair is established by the Coast Guard 
after a secondary station has been installed and the chain is ready to become oper- 
ational. The process to obtain the CSTD begins when the Coast Guard calibrates 
a cesium beam frequency standard, with an optional time standard, at the United 
States Naval Observatory (USNO) or the National Institute of Standards and Tech- 
nology (NIST). The cesium time standard is operated for about a month, tuning it so 
that it will match the master clock, or universal time constant (UTC), at USNO. The 
Coast Guard also tries to minimize the drift of the cesium time standard. [Ref. 9] 

Once the cesium time standard is calibrated, the Coast Guard installs it in a 
master station to use as a time standard. The cesium time standard is compared to 
the time of transmission (TOT) of the master station. The TOT is adjusted to a pre- 
calculated value that is based on the group repetition interval (GRI) and the baseline 
length. Once the TOT is determined, the cesium time standard is brought to each of 
the secondary transmitter sites and connected similarly to that at the master station. 
The secondary TOT is then adjusted to its pre-calculated TOT value. [Ref. 9] 
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During chain calibration, the emission delay (ED) is calculated from the dif- 
ference between the master station TOT and the secondary station TOT; 

EDs = TOTs - TOTm (3.2) 

where the subscripts S and M denote secondary and master, respectively. Once the 
emission delay is within specification, a correction to the CSTD is calculated: 

COR = EDp — EDs (3.3) 

where EDp is the published ED for the baseline. The correction is then added to the 
presently assigned CSTD to calculate the value for the new CSTD. This value then 
becomes the CSTD for a particular master-secondary baseline. CSTDs only exist for 
A1 monitor sites. The A2 monitor sites do not have a CSTD. [Ref. 9] 

2. Phase Time Difference 

The phase time difference (TD) is the difference between the TOA of a master 
station’s signal and the TOA of a secondary station’s signal from the same group repi- 
tition interval (GRI). Continuous measurement of the TD of each respective master- 
secondary pair is made at a monitor location, or locations, within the defined service 
area. This TD is maintained at the CSTD by the insertion of LPAs. In general, the 
hourly average of TD is held to within ±50 ns of the CSTD. 

E. CALOC PROJECT 
1. Background 

The Loran-C control problem is one of controlling a TD, specifically, the differ- 
ence between the times of arrival (TOA) of a master station and a secondary station 
signal. In a study conducted in the mid-1970s, data was collected at the USCG 
Electronics Engineering Center (EECEN) over three unique propogation paths. The 
path from EECEN to the M station was nearly all overland, coastal direct path; 
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Y station was entirely over water; and Z station contained mid-western plain and 
the Appalachian mountains through Pennsylvania. In order to study the nature of 
the short term time difference (TD) fluctuations, as well as model them, TOAs were 
measured for M-Y and M-Z. This process enabled the TD to be measured more ac- 
curately since the TOA for each station was compared to its local cesium oscillator. 
[Ref. 4] 

The TOA records used to develop the model that became the existing con- 
troller were generated from late October to December 1974. Correlation analysis 
performed on the data sets determined that a slowly-varying, long-term component 
was present in the data. After determining that the source of the long-term correla- 
tion was due to a long-term TOA fluctuation, it was approximated as a linear trend 
and removed using linear regression. The resulting autocorrelation with the linear 
trend removed exhibited a more expected behavior. The significant result of the TOA 
correlation analysis was that the TOA fluctuation is composed of two components: a 
short-term component with correlation of less than two hours and an rms value on 
the order of 20-30 ns, and a long-term component that appears as a linear trend over 
periods of 12-24 hours with a slope on the order of 10-20 ns/hour. [Ref. 4] 

Two experiments were conducted to attempt to isolate the possible sources of 
the short-term TOA fluctuations. Due to the propagation effects being distributed 
along lengthy baselines, it was expected that the fluctuations would be due to trans- 
mitted signal synchronization. One source of fluctuation was the local cesium os- 
cillator, which maintains real-time synchronization in each LORSTA. The second 
major source of fluctuation was the cycle-compensation servo loop, which measures 
the pulse crossover timing against a reference derived from the oscillator. The source 
of the fluctuation was the 12-second servo loop time constant, which provided the 
fastest transient recovery and also introduced the greatest amount of noise due to the 
random-phase- aver aging measurement technique. [Ref. 4] 
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An attempt was made to quantify the fluctuations. It was determined that 
the primary source of TOA fluctuations were the random phase noise of the local 
cesium oscillators with a magnitude on the order of 7-14 ns. Additionally, the cycle- 
compensation servo loop did not contribute a statistically significant amount to TOA 
fluctuation. However, the servo loop and other unknown effects did make up some 
portion of the observed fluctuations. [Ref. 4] 

2. Models 

From the analysis conducted in the EECEN study as outlined above, it was 

determined that the Loran-C signal TOAs and, correspondingly, the TDs displayed 

fluctuations having two distinct characteristics: a short-term fluctuation with minimal 

correlation beyond the five minute sampling period used during the study; and a 

long-term fluctuation which appears as a linear trend over a 12-24 hour period. The 

data was insufficient to determine the source of the long-term fluctuation. This was 

deemed unimportant since control of the TDs would be over a much shorter period. 

Two models were explored to provide TD control: a minimum linear model and a 

more complex curve-fit and smooth approach. [Ref. 4] 
a. Linear Model 

The minimum linear model, which accounted for only the short term 
TOA fluctuation, did not make use of any smoothing or averaging techniques. It sim- 
ply used the current sample value as the present best estimate for use in determining 
the correct control signal. The model was inadequate for several reasons: it did not 
recognize the contribution of the long-term TOA fluctuation nor the rapid changes in 
received signal-to-noise ratio (SNR) that may occur due to heavy weather activity; it 
was not desired to generate a control signal in quanta of 20 ns every sample period; 
and it would not adequately control the TD for received SNR conditions significantly 
worse than the observations recorded for this study. [Ref. 4] 

The limitations of a minimum linear model meant that the control al- 
gorithm would not make a major reduction of the short-term TOA fluctuation since 
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that would lead to an undesireable return to fully synchronized operation. Thus, a 
short-term TD fluctuation of 20-30 ns was the best that could be attained with the 
existing equipment and a new algorithm. In order to include the long-term fluctua- 
tions, a model consisting of an integrator driven by white noise with a pole on the 
unit circle at z=l was used. [Ref. 4] 

This estimator of TOA fluctuation was a simple form of Kalman fil- 
ter, as depicted in Figure 6. The value for K was determined empirically from the 
recorded TOA values. However, the value of K could have been determined for ei- 
ther a transient or a steady-state condition given the variance of the process driving 
the integrator and the variance of the received atmospheric noise. These parameters 
represent the time- varying nature of the TD control problem, and therefore a unique 
value could not be determined. [Ref. 4] 




Figure 6. Estimator for Single Integrator Source Model. [Ref. 4] 
b. Curve- fit Model 

An alternative approach to the TD control problem explored a form of 
curve-fitting the previous TOA samples with an n*^^ order polynomial. There were 
several parameters considered in selecting the proper curve-fit model: the fit interval, 
the type and degree of curve to fit, and the fitting method. The biggest advantage 
of curve fitting was that it required no a priori knowledge of the TOA fluctuation 
process or the additive received noise, the parameters that were highly variable with 
time and location. [Ref. 4] 
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The choice of a 3'^'* order polynomial fit over several hours was purely 
subjective, based on qualitative observation and the observer’s judgment and expe- 
rience. The quadratic curve-fit of the TOA data provided one additional degree of 
freedom over its generally linear representation. The selected length of time over 
which to smooth the TOA data was 1^ hours, or 18 data points. The 3''^ order poly- 
nomial fit provided an offset, linear trend, and a quadratic curve to smooth the TD 
data. The least squares approach taken was expressed in matrix form as; 
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where S is an array of received TD samples with Sj\! the most recent sample; Z is an 
array of smoothed TD estimates with the most recent estimate; and H is an Ax 3 
third order polynomial matrix. [Ref. 4] 



F. CALOC ALGORITHM 

The comparison of the linear model and the curve-fit model failed to defini- 
tively prove that one was superior to the other. The curve-fit model was selected based 
on the improvement in using smoothed data versus real-time estimates and its supe- 
rior adaptability to changing TD fluctuation and atmospheric noise conditions. The 
superior adaptability is primarily a consequence of the model being non-recursive: any 
transient disturbances beyond the smoothing interval have no effect while a recursive 
filter, like the linear model, has an infinite memory to these transients. [Ref. 4] 
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1. Performance Metric 

The USCG has established control procedures for how the Loran-C system 
will be operated. These control procedures may be summarized as: 

1. Insert LPAs to maintain the current time difference error (TDE) at the control 
station to of the tolerance bound, which is nominally ± 100 ns; 

2. Minimize the number of LPAs; and 

3. Maintain the cumulative TD close to the CSTD.fRef. 4j 



The rationale behind these procedures is related to the accuracy and the stability 
of the Loran-C system. The predictable accuracy^ is enhanced by maintaining the 
current TDE near zero, and the repeatable accuracy^ is enhanced by maintaining 
the cumulative TD average at the CSTD. [Ref. 1] Minimizing the number of LPAs 
will provide a more stable system by not responding immediately to instantaneous 
fluctuations in noise level due to rapidly changing atmospheric conditions or other 
interference. [Ref. 4] 

2. Cost Function 

In order to be able to quantify the accuracy of any algorithm, a control law 
must be defined. The Loran-C control procedures were translated into a cost function: 



j(„) = ~ + 



u 



(3.8) 



/C| ' K} ' (20ns)2 

where u is the value of the individual LPA being evaluated, Nq is the control num- 
ber, and e{u) is the total predicted cumulative error. Qualitatively there is a cost 
associated with allowing the TD to deviate from the control number, Nc, and there 
is also a cost to correct the error with an LPA. Similarly, there is a cost associated 



'Predictable accuracy is the accuracy of a position with respect to the geographic or geodetic 
coordinates of the earth. [Ref. 1] 

^Repeatable accuracy is the accuracy with which a user is able to return to a known position whose 
coordinates have been measured at a previous time with the same navigational system. [Ref. 1] 
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with the cumulative TDE. The third term in the cost function is a penalty for using 
large LPA corrections. [Ref. 4] 

The terms Kg and Ki are parameters used to tune the control functionality. 
Reducing the value of Kg causes any TD deviation from Nc to be costly while reduc- 
ing Ki places the control emphasis on maintaining the cumulative TDE near zero. 
Increasing both Kg and Ki will cause the cost of an LPA to be increased such that few 
LPAs will be recommended until either the current TDE, {Z^ — Nc), or cumulative 
TDE becomes excessively large. [Ref. 4] 

In this chapter, a brief description of CALOC was presented with a focus on the 
task of time difference control. The significant portions of the original CALOC study 
were presented, including the existing algorithm and associated cost function. The 
next chapter will explore the proportional-integral-derivative (PID) controller. The 
development process of the PID controller includes preprocessing the recorded data 
and testing the algorithm. This, along with the results of testing, will be presented 
in the next chapter. 
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IV. PROPORTIONAL-INTEGRAL- 
DERIVATIVE CONTROL ALGORITHM 



In this chapter, an algorithm based on the proportional-integral-derivative 
(PID) controller is developed. The input is actual data from operational Loran-C 
chains and the output values are local phase adjustments (LPAs). These LPAs are 
used to help maintain the station pulse transmission timing in line with the USCG 
established control procedures. The data provided is preprocessed in order to remove 
the influence of the current controller, CALOC, to obtain uncontrolled data. The data 
is filtered to remove the high frequency noise components and then passed through the 
PID controller. The output of the controller is quantized to produce LPAs in quanta 
of 20 ns for pulse timing correction. The results of this controller show improved 
control of both the current time difference error (TDE) and the cumulative TDE over 
the existing control algorithm. 

A. TIME DIFFERENCE ERROR (TDE) DATA SETS 

The data sets used in this research were collected during three distinct times 
of the year. The data from Malone, Florida, was collected during the winter; data 
from Middletown, California, was collected during the late spring; and the data from 
Kodiak, Alaska, was collected during the fall. This gave some indication of seasonal 
weather effects, however, no conclusions could be drawn from the limited data sets 
and the fact that each chain’s data set included at most two days’ worth of data. 

1. Data Set Locations 

There are eight Loran-C chains in North America that are operated by the 
USCG; North Pacific (NORPAC); Gulf of Alaska (GOA); U.S. West Coast (USWC); 
North Central U.S. (NOCUS); South Central U.S. (SOCUS); Great Lakes (GLKS); 
South-East U.S. (SEUS); and North-East U.S. (NEUS). The data used in this research 
encompasses five of the eight chains, though the data for NORPAC was recorded on 
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one day for only 15 hours. In the following discussion, the chain locations will be 

characterized, and the required preprocessing will be described. 

a. Malone, Florida 

Data collected from Malone, Florida was recorded from two chains: 
South-East United States (SEUS) and South-Central United States (SOCUS). It was 
recorded over a three-day period during 18-20 January 1997, though only the data 
from 18 and 20 January were used due to corruption of the 19 January data. The 
SEUS chain has a master station and four secondary stations with locations listed in 
Table I. The M-X baseline is primarily over water across the Gulf of Mexico while 
the other three baselines (M-W, M-Y and M-Z) are over coastal land masses. The 
estimated ground wave coverage area is depicted as the dashed line in Figure 7. 



Station 


Location 


M 


Malone, Florida 


W 


Orangeville, Louisiana 


X 


Raymondville, Texas 


Y 


Jupiter, Florida 


Z 


Carolina Beach, North Carolina 



Table I. South-East United States Loran-C Chain Transmitter Locations. [Ref. 1] 

The SOCUS chain has a master station and five secondary stations 
with locations listed in Table II. All five master-secondary baselines are over land, 
generally in the central plains area of the United States. The estimated groundwave 

coverage area is depicted as the dashed line in Figure 8. 

b. Middletown, California 

The data collected from Middletown, California was recorded on 
4-5 May 1997 from two chains: United States West Coast (USWC) and North-Central 
United States (NOCUS). The USWC chain has a master station and three secondary 
stations with locations listed in Table III. The master-secondary baselines are all 
over land in the Rocky Mountain range of the western United States. The estimated 
groundwave coverage area is depicted as the dashed line in Figure 9. 
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Figure 7. South-East United States Loran-C Chain. Transmitter Stations: M - 
Malone, FL; W - Grangeville, LA; X - Raymondville, TX; Y - Jupiter, FL; Z - 
Carolina Beach, NC. [Ref. 1] 



1200 w 100* W 80* W 




Figure 8. South Central United States Loran-C Chain. Transmitter Stations: M - 
Boise City, OK; V - Gillette, WY; W - Searchlight, NV; X - Las Cruces, NM; Y - 
Raymondville, TX; Z - Grangeville, LA. [Ref. 1] 
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Station 


Location 


M 


Boise City, Oklahoma 


V 


Gillette, Wyoming 


W 


Searchlight, Nevada 


X 


Las Cruces, New Mexico 


Y 


Raymondville, Texas 


Z 


Orangeville, Louisiana 



Table II. South Central United States Loran-C Chain Transmitter Locations. [Ref. 1] 



Station 


Location 


M 


Fallon, Nevada 


W 


George, Washington 


X 


Middletown, California 


Y 


Searchlight, Nevada 



Table III. United States West Coast Loran-C Chain Transmitter Locations. [Ref. 1] 

The NOCUS chain also has a master station and three secondary sta- 
tions with locations listed in Table IV. The master-secondary baselines in this chain 
are located in the northern plains region of the United States and southern Canada. 

The estimated groundwave coverage area is depicted as the dashed line in Figure 10. 
c. Kodiak, Alaska 

The data collected from Kodiak, Alaska was recorded on 11 Septem- 
ber 1996 for a period of 15 hours from Loran-C Station Zulu of the North Pacific 



Station 


Location 


M 


Havre, Montana 


W 


Baudette, Minnesota 


X 


Gillette, Wyoming 


Y 


Williams Lake, Canada 



Table IV. North Central United States Loran-C Chain Transmitter Locations. [Ref. l] 
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Figure 9. United States West Coast Loran-C Chain. Transmitter Stations: M - 
Fallon, NV; W - George, WA; X - Middletown, CA; Y - Searchlight, NV. [Ref. 1] 




Figure 10. North Central United States Loran-C Chain. Transmitter Stations: M 
- Havre, MT; W - Baudette, MN; X - Gillette, WY; Y - Williams Lake, Canada. 
[Ref. 1] 
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(NORPAC) Chain. The NORPAC chain has a master station and three secondary 
stations with locations listed in Table V. The master-secondary baselines are primar- 
ily over water in the Bering Straits, though the M-Z baseline passes over the Alaskan 
peninsula. The estimated groundwave coverage area is depicted as the dashed line in 
Figure 11. 



Station 


Location 


M 


St. Paul, Alaska 


X 


Attu, Alaska 


Y 


Port Clarence, Alaska 


Z 


Kodiak, Alaska 



Table V. North Pacific Loran-C Chain Transmitter Locations. [Ref. 1] 
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Figure 11. North Pacific Loran-C Chain. Transmitter Stations: M - St. Paul, AK; X 
- Attu, AK; Y - Port Clarence, AK; Z - Kodiak, AK. [Ref. l] 



2. Limitations 

There were several limitations in the data sets used in this work. The most 
significant of these was that the data sets did not contain 24 hours worth of data 
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points. The initial data sets from Malone, Florida were missing anywhere from 65 
minutes to as much as 164 minutes of data points. This precluded an exact corre- 
lation with the CALOC-generated strip charts and forced a different approach for 
comparison. After isolating and correcting the cause of the missing data points, the 
data sets from Middletown, California, were only missing between 7 and 17 minutes 
worth of data. Due to this limitation, each data set was used to generate strip charts 
that included only the data points. 

A second limitation was that the data from each location (Malone, Florida, 
Middletown, California, and Kodiak, Alaska) covered no more than two days during 
the same time of the year. No comparison could be made to determine seasonal 
trends nor to verify that the data received was normal for each time of year. In 
order to perform a rigorous treatment of the data and to implement a controller, a 
significant amount of data would need to be collected during different seasons and 
under varying atmospheric conditions. Since this was not the case, additional testing 
on actual equipment under varying conditions is required prior to implementation of 
the proposed algorithm. 

A third limitation was that the data sets included the effects of CALOC. Due 
to Loran-C’s significant role in coastal navigation, there was no acceptable method 
to collect uncorrected data from an operational Loran-C chain; without CALOC’s 
control input, the Loran-C chain would transmit an inaccurate signal during the time 
of data colletion. This required deriving a method for removing as much of CALOC’s 
influence as possible so that the new contoller would have raw data that only included 
the inherent noise due to the propagation path, man-made interference, atmospheric 
noise, etc. 

3. Data Preprocessing 

There was little that could be done to mitigate the limitations of incomplete 
data sets and limited seasonal variation; however, it was believed that the effects of 
CALOC on the data set could be removed. The data sets, once plotted, revealed a 
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significant delay between the time CALOC ordered an LPA and the time that its full 
effect was detected at the monitoring station. Thus, the first step in removing the 
effect of CALOC on the data was to determine the average time of this processing 
delay. Data preprocessing for each data set was accomplished using the Matlab code 
contained in Section 5 of Appendix D. After loading the data set, the preprocessing 
involves manually determining the timing of LPAs within the data set and then 

removing them using a time-phased correction described below. 

a. Processing Delay 

Based on the strip charts generated by CALOC that were provided 
with the data sets, it appeared that the time delay experienced between CALOC 
ordering a correction and its effect being detected at the monitor station was relatively 
constant. As a first order approximation, it was assumed that this delay was less than 
the 7.5 minute sampling interval presently used in CALOC. However, upon further 
investigation, it was found that the average time delay was on the order of 56 data 
points, or over nine minutes. Identification of the source of this time delay is beyond 
the scope of this work. 

To determine the magnitude of the time delay, a moving average over 
45 data points, equating to CALOC’s 7.5 minutes average interval, was taken of each 
data set and plotted. The time delay was displayed in the graph by a vertical change 
in the time difference error (TDE) as may be seen in Figure 12. The time delay was 
determined by counting the number of 10-second samples between the indication in 
the data of an LPA and the indication that the pulse transmitter had shifted to the 

new transmission time. 

b. Corrected Data 

Based on 146 LPAs in the nine data sets, the average time delay was 
computed to be 556.4 seconds. In order to apply the correction to remove the LPAs 
due to CALOC, the time delay was rounded up to 560 seconds, which equated to 
56 10-second data points. The LPAs were then removed by applying the inverse 
of the LPA uniformly over the 56 data points. The same data described above, 
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Figure 12. South-East United States Loran-C Station Yankee Time Difference Error 
recorded on 18 January 1997. 




Figure 13. South-East United States Loran-C Station Yankee Time Difference Er- 
ror recorded on 18 January 1997 with CALOC Local Phase Adjustments (LPAs) 
Removed. 
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after removing the LPAs, may be seen in Figure 13. The time delay displayed no 
dependence on the magnitude of the ordered LPA since it was consistent for both a 
±20 ns LPA or a ±40 ns LPA. This indicated that the delay was due solely to the 
time required to apply the LPA to the pulse transmitter. 



B. FINITE IMPULSE RESPONSE (FIR) FILTER 

The recorded Loran-C data contains high frequency noise components that 
need to be removed. These are due primarily to noise sources in the wireless medium, 
but also include noise from the transmitter and receiver. Additional noise is injected 
in the data during the preprocessing. To remove these high frequency noise sources, 
a digital low-pass filter is required. There are two basic types of digital filters: finite 
impulse response (FIR) and infinite impulse response (IIR). An FIR filter was selected 
primarily due to its linear phase characteristic. Other benefits of the FIR filter are 
that it is simpler to design and, when windowed, is of finite duration. 

There are several methods of FIR filter design that may be used, including 
frequency sampling and windowing. The windowing method truncates an infinite 
duration unit sample response, hd{n), to length M — 1 by multiplying hd{n) by a 
rectangular window of length M, defined as: 



w{n) = 



1, n = 0, 1, . . . , M — 1 



(4.1) 



0, otherwise. 

A symmetric lowpass linear-phase FIR filter is described in the frequency domain as: 



HdH = ( 






0 < Ittll < UJc 



0, otherwise. 

A windowed FIR filter of length M has a unit sample response 



(4.2) 



h{n) = 



sino)c(n — 



0 < n < M — l,n ^ 



M -1 



(4.3) 
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where Wc is the cut-off frequency and M is the filter length. Using a rectangular 
window has one significant drawback; it causes relatively large sidelobes. To elimi- 
nate the sidelobes, a window function that gradually decays toward zero is desirable. 
[Ref. 10] 

There are several windowing options available. Here we use the Hamming 
window, which is given as: 

/ , . 2nn , 

tu(n) = 0.54 — 0.46 cos(-^^^^ — j-), 0<n<M — 1 (4.4) 

where Mis the window length or filter length. The Hamming window has significantly 
lower sidelobes compared to the rectangular window. However, for the same length 
filter, the width of the main lobe is wider which leads to a wider transition region in 
the FIR filter response. [Ref. 10] The filter used in this application is an FIR filter 
with a Hamming window of length 32 and a cut-off frequency of 7 t/ 45. The cut-off 
frequency was selected based on the removal of all high frequency components beyond 
the 45 data points used in each average. 

The FIR filter is generated using the function firham.m (see Section 2, Ap- 
pendix D). The inputs to the function are the window length, M, and the cut-off 
frequency, cvc', the outputs are the filter coefficients. For M = 32, Uc = tt/ 45, and 
using a Hamming window, the coefficients of the FIR filter are listed in Table VI. 



h(0) = h(31) 


0.0040 


h(l) = h(30) 


0.0046 


h(2) = h(29) 


0.0062 


h(3) = h(28) 


0.0088 


h(4) = h(27) 


0.0123 


h(5) = h(26) 


0.0166 


h(6) = h(25) 


0.0217 


h(7) = h(24) 


0.0272 



h(8) = h(23) 


0.0330 


h(9) = h(22) 


0.0389 


h(10) = h(21) 


0.0445 


h(ll) = h(20) 


0.0496 


h(12) = h(19) 


0.0540 


h(13) = h(18) 


0.0575 


h(14) = h(17) 


0.0600 


h(15) = h(16) 


0.0612 



Table VI. FIR Filter Coefficients. 
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The response of the FIR filter prior to windowing with the Hamming window is 
depicted in Figure 14. The peak of the first sidelobe is approximately -20 dB below the 
passband peak, and the phase response is linear over the passband. After windowing, 
the filter response, shown in Figure 15, has a peak sidelobe level of approximately 
-45 dB while maintaining the linear phase. As mentioned above, the transition band 
has increased when using the Hamming window compared to the rectangular window. 





Figure 14. Magnitude and Phase Response for Finite Impulse Response (FIR) Fil- 
ter (Length = 32, Cut-off Frequency = tt/ 45) with no Windowing (i.e., rectangular 
window). 



C. PROPORTIONAL-INTEGRAL-DERIVATIVE (PID) 
CONTROLLER 

This section presents a PID controller. The process began with exploring a 
model of the existing control system and implementing that model. Based on the 
model, suitable weights were derived for the controller, and then the algorithm was 
tested using all the data sets. 
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Figure 15. Magnitude and Phase Response for Finite Impulse Response (FIR) Filter 
(Length = 32, Cut-off Frequency = tt/ 45) with Hamming Window. 

1. Theory 

The PID controller is widely used in industrial processes to control systems 
that cannot be modeled as lumped and linear time-invariant. It is essentially three 
controllers in one: a proportional controller, an integral controller and a derivative 
controller. [Ref. 11] 

The transfer function of the proportional controller is a gain, kp. For a given 
controller input, e{t), the output will be: 

u{t) = kpe{t). (4.5) 

For a stable plant, the feedback system will remain stable for a range of kp. However, 
as kp increases, the unit-step response may become faster and the feedback system 
becomes unstable. One significant limitation to the proportional controller is that for 
the same unit-step reference input, the steady-state plant output will be different for 
different values of kp requiring the plant set point to be manually reset. [Ref. 11] 
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For the same input, the integral controller output is: 



u{t) = ki [ e{T)d,T 
Jo 



(4.6) 



where A;, is the integral constant. This leads to a transfer function for the integral 
controller of For a stable feedback system, the steady-state error due to any step- 
reference input is zero for all ki. This means that there is no need to manually reset 
the set point when k{ is changed. For this reason, integral control is also called reset 
control with ki being the reset rate. Though the requirement to reset the set point has 
been eliminated, the integral controller makes stabilizing the feedback system more 
difficult. [Ref. 11] 

The derivative controller maintains control of a system by using the rate of 
change of the error signal as its input. The output of the derivative controller is: 



u{t) = kd 



de{t) 

dt 



(4.7) 



where e{t) is the input to the controller and k^ is the derivative constant. The 
derivative controller by itself is rarely used in feedback control systems because if the 
error signal changes slowly or is constant, the controller output would be minimal or 
zero allowing a potentially large error to remain. [Ref. 11] 

A PID controller simultaneously takes advantage of all three controllers de- 
scribed above. The three parameters of the controller, kp, ki and kd, are tuned to 
achieve the desired system response. One method of tuning is by trial-and-error, ini- 
tializing each parameter to some constant value and varying all three parameters to 
achieve a satisfactory system output. Exact values of these parameters could not be 
obtained due to lack of an accurate model for the entire Loran-C control subsystem. 
Because of this, a trial-and-error approach was used to obtain the solution described 
in this chapter. [Ref. 11] 
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2. Model 

The most important step in the process of developing a control algorithm is 
constructing an accurate model. The system model is a mathematical representation 
of the physical system. It may be used to evaluate the performance of the physical sys- 
tem to different input stimuli when the system is unavailable for operational testing, 
as is the case for the Loran-C system. The model developed for the existing control 
system, shown in Figure 16, contains four major sections, each modelled as separate 
transfer functions: CALOC, the transmitter, the medium and the monitoring station. 




Figure 16. Loran-C Control System Model. 

Knowledge of the transfer functions for the transmitter, medium and the mon- 
itoring station is not required based on the method of data preprocessing discussed 
earlier. The preprocessing only removed the effect of CALOC while leaving intact all 
noise components and the dynamic behavior of the transmitter, the medium, and the 
monitoring station equipment. This significantly simplified the process of modelling 
the new control system. All that was required was to substitute the new control 
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transfer function into the block that contained the CALOC transfer function. The 
new model is depicted in Figure 17. 




Figure 17. Representative Model of the Loran-C Control System with PID Controller. 

Some modifications to the model were required to support simulation of the 
new control system. A method to input local phase adjustments (LPAs) in the data 
was developed to simulate the functionality of the new control algorithm within the 
Loran-C control system. As depicted in Figure 18, the simlation model sums the 
generated LPA with the future time difference errors (TDEs) from the recorded data. 
The LPA process was modelled as an integrator and time delay with a transfer func- 
tion of: 



Gl(z) = (4.8) 

The LPA is a timing adjustment to the Loran-C station pulse transmitter for all 
future transmissions which leads to modelling the LPA function as an integrator. 
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Figure 18. Simulation Model of the Loran-C Control System. 

The single time delay was based on the fact that the effect of the LPA was detected 
at the monitoring station in about one 7.5-minute control interval. 

The algorithm for the PID controller is the classical model shown in Figure 19. 
The input to the controller is the TDE at time k and the output is the continuous 
time LPA recommendation for time k. The proportional gain, A'p, is applied to the 
TDE directly. The integral of the TDE is the cumulative TDE and is controlled using 
the integral gain, Ki. The error rate of change is calculated by taking the difference 
between e(/c) and e{k - 1). The derivative gain, Kd, is applied to this difference. 
Summing all three contributions of the PID controller yields the LPA output that is 
then provided as input to the quantizer, which yields discretized values, u{k). 

3. Quantization 

Quantization is most commonly used in analog-to-digital (A/D) conversion to 
convert a continuous time signal into a set of discrete values. Here, quantization is 
used to transform the LPA output of the controller into quanta of 20 ns, over a range 
of ±180 ns, for use in controlling pulse transmission timing. There are two basic types 
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Figure 19. PID Controller without Quantization. 

of quantizers: mid-rise and mid-tread. The mid-tread quantizer is generally preferred 
over the mid-rise because the quantization noise has zero mean. In this application, 
the mid-tread is used as it did not introduce a bias to the signal. [Ref. 10] 

Quantization does not come without some cost. The process of quanitization 
introduces an error that cannot be removed. Figure 20 depicts the uniform, also 
known as linear, quantizer designed for this controller. The step size is a uniform 
value of 20 ns, which yields a maximum quantization error of 10 ns, or half the step 
size, throughout the dynamic range of the quantizer. The quantization error for each 
sample may be modeled as a zero-mean, uniform random variable over the interval 
±10 ns. 

D. PID RESULTS 

The PID control algorithm was implemented using the function pid.m (see 
Section 3, Appendix D). The inputs to the function are the data set name, initial 
condition for the cumulative TDE, and title for the plot. The plot title is generated 
by the function parse.title.m (see Section 7, Appendix D). 
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Figure 20. PID Controller Output Quantizer for Generating Local Phase Adjust- 
ments. 

The PID control algorithm proved to be a definite improvement over the curve- 
fit algorithm in the existing Loran-C controller. In every case, both the current time 
difference error (TDE) and the cumulative TDE are maintained in closer control. In 
some cases, the PID controller was even able to accomplish this with fewer local phase 
adjustments (LPAs) than CALOC. However, in eight cases out of 29 there are a large 
number of LPAs being ordered which is counter to one of the USCG control policies: 
to minimize the number of LPAs. This excessive chatter of the controller could be 
minimized with further operational testing using a greater number of data sets. 

The results are presented for two Loran-C stations: SEUS Station Yankee and 
SEUS Station Zulu. The results of the other 27 data sets are included in Appendix B. 
SEUS Station Yankee cumulative TDE was initialized at -31 ns based on the actual 
value found on the CALOC strip chart. As is shown in Figure 21, CALOC maintained 
the current TDE within a range of -37 ns to 28 ns while the cumulative TDE has a 
range of ±50 ns. The PID controller showed a 14% improvement in dynamic range 
in both current TDE, with a range of -30 ns to 26 ns, and cumulative TDE, with a 
range of -40 ns to 46 ns. 
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The results for SEUS Station Zulu were even better. The cumulative TDE 
was initialized at -77 ns. As seen in Figure 22, CALOC was able to maintain the 
current TDE between -38 ns and 32 ns and the cumulative TDE between -71 ns and 
38 ns. The PID controller was able to provide a significantly smaller dynamic range 
of cumulative TDE for Station Zulu. For the PID controller, the current TDE had 
a range of ±31 ns, an 11% improvement, and the cumulative TDE was maintained 
between -53 ns and 30 ns for a 24% improvement. 

The results of the PID controller on the 29 data sets indicate that it is a 
viable option to CALOC. There was improved control of the TDE in every case 
though some data sets required a greater number of LPAs than the same data under 
CALOC. However, in the case of the data from the United States West Coast (USWC) 
Chain and the North Central United States (NOCUS) Chain, the chains are rarely in 
automatic control, which allows CALOC to function fully with no human intervention. 
Thus, the direct comparison with CALOC on those chains is not as conclusive. 

Based on the overall results and the knowledge of the ever-changing medium in 
which Loran-C operates, the control algorithm needs to be data-dependent to better 
respond to this high noise environment. A PID controller could provide the needed 
adaptability; however, there are other controllers that are better suited. One of these 
is the Kalman filter, which is developed in the next chapter. 
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TDE (ns) TDE (ns) 




(a) CALOC 




(b) PID Controller 

Figure 21. Strip Chart for South East United States (SEUS) Loran-C Chain Station 
Yankee recorded on 18 January 1997; o Linear Phase Adjustment (LPA) (magnitude 
times 2); + Time Difference Error (TDE); — Cumulative TDE. 
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(a) CALOC 




(b) PID Controller 



Figure 22. Strip Chart for South East United States (SEUS) Loran-C Chain Station 
Zulu recorded on 18 January 1997: o Linear Phase Adjustment (LPA) (magnitude 
times 2); + Time Difference Error (TDE); — Cumulative TDE. 
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V. DATA-DEPENDENT CONTROL 

ALGORITHM 

The PID controller was an initial step in determining whether significant im- 
provement over CALOC was possible. As was proved in the preceding chapter, there 
is significant room for improvement. This chapter will demonstrate the effects of a 
more robust controller, the Kalman filter. One of the advantages of the Kalman fil- 
ter is its ability to adapt to changes in the operational environment. The discussion 
presented here is introductory, and the results reported are preliminary. 

A. DATA-DEPENDENT SOLUTION 

The primary reason for developing a data-dependent solution for the Loran-C 
control problem is the operating environment of the signal. The signal environment 
is a wireless media that is subject to these sources of interference, both natural and 
man-made. These sources of interference include atmospheric changes, terrain and 
electromagnetic intereference. 

1. Atmospheric Effects 

The sources of atmospheric effect on the Loran-C signal are primarily due to 
temperature changes and seasonal weather effects. The effect of temperature changes 
on signal propagation has no reasonable theoretical explanation other than its effect 
on the refractive index of the media. Specifically, the Loran-C signal exhibits a 
diurnal variation that can be associated with temperature changes. [Ref. 12] One 
source of these temperature effects is a temperature inversion and associated humidity 
change which may produce highly variable refractivity profiles. This process produces 
an evaporation duct in which radio waves may propagate well beyond the horizon. 
Overland propagation refraction effects may be seen in at least two situations; one 
is the evaporation of water on the ground after raining which produces an over-land 
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evaporation duct, and the second is the radiational cooling over deserts and snow- 
covered terrain on clear nights which yields anomalous propagation. [Ref. 13] 

The seasonal weather effects on the Loran-C signal are caused by the passage of 
cold and warm fronts predominate at specific times of the year. The interpretation of 
the correlation between the time difference and total refractivity is that the passage of 
cold and warm fronts produce transient decreases in the propagation time delay. The 
cold front decreases occur after the front has passed while the decreases associated 
with the warm front occur before frontal passage. [Ref. 12] 

The conclusion that may be drawn from these effects is that the Loran-C signal 
is being influenced by changes in the refractive index well above the earth’s surface. 
A more complete understanding of these atmospheric effects is required in order to 
make predictions. In the case of the Loran-C system, accurate predictions correspond 
directly with greater position accuracies. [Ref. 12] 

2. Terrain 

The signal effects of terrain on the Loran-C signal are due to diffraction phe- 
nomena, which are exhibited in both diffuse scattering and coherent scattering. These 
phenomena may be broken down into two categories. The first is specific irregularities 
in the terrain, such as mountain ridges, which are often modeled as knife-edges or, 
in some cases, cylinders. Secondly, diffraction is produced by the curvature of the 
smooth spherical surface of the earth which applies to propagation over the ocean or in 
the plains. At low frequencies, diffraction by hills or ridges alters the field strength in 
a significant region around the line of sight which means that more than one diffract- 
ing feature on the terrain may be contributing to the propagation losses. This is a 
challenging issue that has not yet been successfully modeled for long ray-paths over 
varying terrains. [Ref. 13] 
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3. Electromagnetic Interference 

The sources for electromagnetic interference for the Loran-C system are es- 
sentially of two types: fixed radiators and mobile, or random, radiators. The fixed 
radiators are those that continuously transmit at a constant frequency or transmit at 
consistent times of the day while the mobile radiators randomly transmit their signal 
and maintain no pattern to their transmission. Each master-secondary baseline will 
have a different set of known radiators that may be relatively easily modeled. How- 
ever, determining a suitable model for the mobile radiator may be quite involved and 
location dependent. 

Two fixed radiators that are significant contributors to electromagnetic inter- 
ference are overhead power lines and radio towers. The power lines are low power 
signals that exist within the Loran-C frequency band. This has become of particular 
interest since the expansion of the Loran-C system across the North American con- 
tinent. The signal amplitudes vary as a function of the square of the distance from 
the power line, so their effect is only critical near to the power line, on the order 
of one mile. However, because so many power lines exist, using notch filtering to 
remove the effects is not feasible. Radio towers along the ray-path also contribute 
to the electromagnetic interference seen by the Loran-C signal. One example of this 
is the frequency shift keying found at naval communications stations. These, how- 
ever, unlike the power lines, may be divided using superposition and treated as either 
deterministic or random noise when considering their effects on the time of arrival. 
[Ref. 14] 

One example of the mobile, or random, radiator that causes electromagnetic 
interference with the Loran-C signal is the cathode ray tube (CRT). The raster scan 
displays on computers and test instruments cause interference due to the harmonics 
of the horizontal scan frequency. These frequencies are extremely stable and may be 
isolated once the device is identified. [Ref. 14] There are other sources of random 
electromagnetic interference that would also need to be accounted for. One prominent 
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example is the ham radio operator that may be infrequently transmitting near the 
operating frequency of Loran-C. The mobile radiator presents a unique challenge in 
that their signal, in the aggregate, is a non-deterministic, random process. 

B. KALMAN FILTER 

Kalman filtering provides several advantages over the PID algorithm in solving 
the Loran-C control problem. Its mathematical formulation is described in state-space 
form, and its solution is computed recursively. Thus, each updated state estimate is 
computed only from the previous estimate and the new input data which requires only 
that the previous state estimate be stored. The Kalman filter is also more efficient 
than computing the next state estimate using the entire past observed data at each 
step in the process. [Ref. 15] 

1. Model 

In the general case, the Kalman filtering problem is stated in terms of two 
vectors of random variables: an M-dimensional parameter x(A:) which denotes the 
state of a discrete-time, linear, dynamic system and y(A;) which denotes the observed 
data of the system. The system model may then be fully described by two equations: 
a process equation and a measurement equation. [Ref. 15] 

The process equation is: 

x{k + 1) = A ■ x{k) + 'v(k) (5.1) 

where A is an Mx M state transition matrix that relates the system state at times 
k and k + 1. The Mxl vector v{k) represents process noise and is modeled as a 
zero- mean, white-noise process, ideally Gaussian. [Ref. 15] 

The measurement equation, which describes the observation vector, is: 

z{k) = C ■ x(k) -I- Aw(A;) (5.2) 
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where C is an iVx M measurement matrix. The iVxl vector Aw(k) represents mea- 
surement, or sensor, noise and is modeled as another zero-mean, white-noise process, 
ideally Gaussian. [Ref. 15] 

The Kalman filter attempts to use the past observed data values, consisting of 
y(0), . . . , y(A:), and the model to predict the minimum mean-square estimate of the 
new state, x{k + 1). Thus, the Kalman filter is attempting to solve the equation: 

x{k + l) = A- x{k) + L{k) ■ [z{k) - C ■ x(k)j (5.3) 

where z(k) - C • x(A:) is the correction to be applied and L(A:) is determined in order 
to minimize S || x(k) |p where the error, x(k) = x(k) - x(A:). 

In tailoring the general model to the specific case of the Loran-C control prob- 
lem, depicted in Figure 23, several adaptations were made. In Equation 5.1, the first 
term becomes x(k) -j- A<F(A:) yielding a new process equation of 

x(k + 1) = x(k) -f A$(fc) -f v(k) (5.4) 

where x(k) is the estimate of the current TDE and A4>(A:) is the input LPA. Similarly, 
the first term in Equation 5.2 becomes z(k) -f x(k) yielding a new measurement 
equation of 



z(k -j- 1) = z(k) + x(k) + Aw(k) (5.5) 

where z(k) is the cumulative TDE. 



LPA 





TDE estimate 



Aw{k) 



Figure 23. Control System Process Model. 
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As discussed in Chapter III, the long-term fluctuations of 12-24 hours were of 
lesser consequence since the control of the time difference (TD) would be over a much 
shorter period. A suitable model for the short-term fluctuations in the time difference 
error (TDE) is a random walk with additive white Gaussian noise as depicted in 
Figure 24. The model is characterized by the equation: 

w{k) - Wjn{k) + Aw{k) (5.6) 

where 



w^{k) = Wm{k-l) + v{k-l) (5.7) 

with w{k) and v{k) being zero mean, uncorrelated, white Gaussian noise. 




Figure 24. Random Walk with Additive White Gaussian Noise 



2. Optimal Controller 

Based on the model described above, the Kalman filter estimates the whole 
state as defined by: 

x(*+i) = m + A$(fe) + - m) (5.8) 

p(t+l) = p(i) + Q-gP^ 
z(A; 4- 1) = z(k)+x(k) 



56 



(5.9) 

(5.10) 



where x{k) is the TDE estimate, and z{k) is the estimate of the cumulative TDE. 
The constants Q and R are defined as 

Q = cov(v(A;)) (5.11) 

R = cov(Aw(A;)). (5.12) 

In order to devise an optimal controller, a control law must be defined. As 
discussed in the previous chapter, the USCG has established Loran-C control proce- 
dures concerning the allowed range of time difference and minimizing the number of 
LPAs. The metric that is used to measure the performance of the controller is a cost 
function, defined as: 



J(^) = + a2\z{k)f + |A^(A;)p} (5.13) 

k=0 

where x{k) is the estimate of the current TDE, z{k) is the estimate of the cumulative 
TDE, and A^{k) is the LPA. The cost function is similar to the one defined for 
CALOC in that it applies a cost to a large TDE, a large cumulative TDE, and a large 
LPA. Control of the process is defined by 



A$(A;) = -L(A:) 



x{k) 

z{k) 



(5.14) 



where L(A:) is determined recursively using the Matlab function dlqr, which solves 
the discrete-time linear-quadratic regulator problem and the associated Riccati equa- 
tions. 



C. RESULTS 

The Kalman filter was implemented using the function kalman.m (see Section 
4, Appendix D). The inputs to the function are the data set name, initial condition 
for the cumulative TDE, and title for the plot. The plot title is generated by the 
function parse.title.m (see Section 7, Appendix D). 
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As in the case of the PID controller, the Kalman filter improved the control of 
the time difference error (TDE) over the existing Loran-C controller. It was also an 
improvement over the PID controller and, in most cases, inserted fewer local phase 
adjustments (LPAs). The results detailed below compare the PID control algorithm 
to the Kalman filter as the PID has already been shown to be superior to the existing 
controller. The same two Loran-C stations used in the preceeding chapter are used 
here to illustrate the improvement of the Kalman filter. The results of the other 27 
data sets are included in Appendix C. 

The initial cumulative TDE values for SEUS Stations Yankee and Zulu are 
the same as for the PID controller: -31 ns and -77 ns, respectively. As may be seen in 
Figure 25, the Kalman filter essentially maintained the same range for current TDE 
while improving control of the cumulative TDE. The range for the current TDE for 
the Kalman filter was -30 ns to 29 ns as compared to a range of -30 ns to 26 ns for the 
PID controller. However, the Kalman filter showed a 31% improvement in cumulative 
TDE with a range of -26 ns to 33 ns compared to the PID controller with a range of 
-40 ns to 46 ns. 

The improvement of the results for SEUS Station Zulu were more pronounced. 
Similar to Station Yankee, Station Zulu maintained the current TDE within a similar 
range: -31 ns to 36 ns as compared to ±31 ns for the PID controller. As Figure 26 
illustrates, the range of cumulative TDE for the Kalman filter was -33 ns to 32 ns as 
compared to -53 ns to 30 ns for the PID controller, a 34% improvement. 
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TDE (ns) TDE (ns) 




(a) PID Controller 




(b) Kalman Filter 

Figure 25. Strip Chart for South East United States (SEUS) Loran-C Chain Station 
Yankee recorded on 18 January 1997: o Linear Phase Adjustment (LPA) (magnitude 
times 2); -f- Time Difference Error (TDE); — Cumulative TDE. 
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(a) PID Controller 




(b) Kalman Filter 



Figure 26. Strip Chart for South East United States (SEUS) Loran-C Chain Station 
Zulu recorded on 18 January 1997: o Linear Phase Adjustment (LPA) (magnitude 
times 2); + Time Difference Error (TDE); — Cumulative TDE. 
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VI. CONCLUSIONS 



Loran has a long history of providing accurate navigational information to 
both mariners and aviators and can have a bright future into the next century as a 
complementary radionavigation system to GPS. To ensure this future of Loran-C, its 
accuracy must be enhanced and operational costs reduced. One method of making this 
improvement is to upgrade the transmitters and receivers in each of the LORSTAs, 
but this would be extremely costly. An alternative is to improve the control algorithm. 
Modernizing the control system for Loran-C requires a control algorithm that not 
only maintains the time difference error of the transmitted pulse near zero, but also 
is responsive to changing atmospheric conditions. To this end, this thesis developed 
a model of the existing system and proposed two control algorithms that provided 
substantial improvement over the existing algorithm: one was the PID controller and 
the other was a Kalman filter. 

A. SIGNIFICANT RESULTS 

Modelling Loran-C’s ever-changing operational environment is extremely diffi- 
cult. It includes several pieces of hardware and a variety of noise sources that all have 
unique characteristics. There currently exists no mathematical model of the Loran-C 
control system that could be used to simulate its functionality. Due to this, only the 
control process of inserting Local Phase Adjustments (LPAs) was modelled leaving 
the remaining system effects in place. This simplified the design effort since actual 
data was available that included these effects. 

However, working with actual system data presented unique challenges. It was 
not possible to have a Loran-C chain provide inaccurate navigational signals simply 
to record information for this thesis. The effects of the current controller, CALOC, 
were included with the data and had to be preprocessed. These effects were removed 
by applying a correction to the data over the average time that elapsed before the 



61 



impact of the LPA was felt. This process provided the raw data for testing the 
proposed control algorithms. 

An understanding of Loran-C, and in particular CALOC, was required in 
order to develop a suitable control algorithm. The study of Loran-C focused on one 
of the four tasks of CALOC, time difference control. The four major components of 
this control effort are CALOC, the LORSTA transmitter, the wireless medium, and 
the monitoring station. The effects of the last three components remained after the 
removal of the CALOC LPAs, allowing for the development of an algorithm without 
having to model each component. 

The PID control algorithm blends the strengths of three controllers into one; a 
proportional controller, an integral controller, and a derivative controller. This algo- 
rithm was developed using a basic model and tuned using a trial-and-error approach 
to provide the best control. The results indicate that this algorithm provides more 
accurate pulse transmission timing than CALOC. Once tuned, the PID controller 
provided substantial improvement over the existing algorithm in all cases, though in 
some cases the frequency of timing corrections was excessive. 

The Kalman filter was explored as an alternative to the PID controller. The 
filter provides the capacity to be adaptive in a changing environment and is much 
more robust than the PID controller. This is an extremely important part of what is 
needed in the future control system for Loran-C in today’s technological environment. 
The preliminary results of the Kalman filter showed an even greater improvement over 
the existing controller than the PID controller. Because the Kalman filter will need 
to be more finely tuned to be completely operational, the transition from CALOC to 
the Kalman filter will need an intermediate step, and the PID controller is the best 
choice for this step. 
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B. RECOMMENDATIONS FOR FURTHER STUDY 

Further study of the Loran-C control system will be valuable in at least three 
areas. First, an accurate mathematical model of the Loran-C control system needs to 
be developed. This will provide the opportunity to easily evaluate possible improve- 
ments to Loran-C while minimizing the operational testing required under normal 
situations. Second, the PID controller needs to be fully operationally tested to de- 
termine what further tuning is required prior to implementation. This will provide 
a suitable bridge between what exists today and the controller of the future. The 
PID controller will remove the necessity for the Loran-C station operators to second- 
guess the system. Finally, the Kalman filter needs to be developed further in order 
to become fully operational. What is presented in this thesis is a proof of concept 
and no attempt was made at optimizing this algorithm. To enable optimization of 
the Kalman filter, a wider variety of data sets from all operational Loran-C chains is 
required. The culmination of this work will be to have a Kalman filter implemented 
that can respond to the Loran-C changing operational environment while providing 
optimal navigational information. 
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APPENDIX A. STRIP CHARTS GENERATED 

BY CALOC DATA 

1. SOUTH EAST UNITED STATES (SEUS) 
a. 18 January 1997 




(a) SEUS Station Yankee 




(b) SEUS Station Zulu 

Figure 27. CALOC Strip Chart for South East United States (SEUS) Loran-C Chain 
recorded on 18 January 1997: o Linear Phase Adjustment (LPA) (magnitude times 
2); + Time Difference Error (TDE); — Cumulative TDE. 
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b, 20 January 1997 




(a) SEUS Station Whiskey 




(b) SEUS Station Xray 
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(c) SEUS Station Yankee 




(d) SEUS Station Zulu 



Figure 28. CALOC Strip Chart for South East United States (SEUS) Loran-C Chain 
recorded on 20 January 1997: o Linear Phase Adjustment (LPA) (magnitude times 
2); + Time Difference Error (TDE); — Cumulative TDE. 
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2. SOUTH CENTRAL UNITED STATES (SOCUS) 
a. 18 January 1997 




(a) SOCUS Station Victor 




(b) SOCUS Station Whiskey 
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(c) SOCUS Station Xray 




(d) SOCUS Station Yankee 




(e) SOCUS Station Zulu 



Figure 29. CALOC Strip Chart for South Central United States (SOCUS) Loran-C 
Chain recorded on 18 January 1997: o Linear Phase Adjustment (LPA) (magnitude 
times 2); + Time Difference Error (TDE); — Cumulative TDE. 
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b. 20 January 1997 
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(a) SOCUS Station Victor 




(b) SOCUS Station Whiskey 
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(c) SOCUS Station Xray 




Time (hours) 

(d) SOCUS Station Yankee 




(e) SOCUS Station Zulu 



Figure 30. CALOC Strip Chart for South Central United States (SOCUS) Loran-C 
Chain recorded on 20 January 1997: o Linear Phase Adjustment (LPA) (magnitude 
times 2); + Time Difference Error (TDE^ — Cumulative TDE. 
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3. NORTH CENTRAL UNITED STATES (NOCUS) 
a. 4 May 1997 




(a) NOCUS Station Whiskey 




(b) NOCUS Station Xray 
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(c) NOCUS Station Yankee 



Figure 31. CALOC Strip Chart for North Central United States (NOCUS) Loran-C 
Chain recorded on 4 May 1997: o Linear Phase Adjustment (LPA) (magnitude times 
2); + Time Difference Error (TDE); — Cumulative TDE. 
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b. 5 May 1997 




(a) NOCUS Station Whiskey 




(b) NOCUS Station Xray 
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(c) NOCUS Station Yankee 



Figure 32. CALOC Strip Chart for North Central United States (NOCUS) Loran-C 
Chain recorded on 5 May 1997; o Linear Phase Adjustment (LPA) (magnitude times 
2); + Time Difference Error (TDE); — Cumulative TDE. 
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4. UNITED STATES WEST COAST (USWC) 
a. 4 May 1997 
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(b) USWC Station Xray 
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(c) USWC Station Yankee 

Figure 33. CALOC Strip Chart for United States West Coast (USWC) Loran-C 
Chain recorded on 4 May 1997: o Linear Phase Adjustment (LPA) (magnitude times 
2); + Time Difference Error (TDE); — Cumulative TDE. 
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b. 5 May 1997 




(a) USWC Station Whiskey 




(b) USWC Station Xray 
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(c) USWC Station Yankee 

Figure 34. CALOC Strip Chart for United States West Coast (USWC) Loran-C 
Chain recorded on 5 May 1997: o Linear Phase Adjustment (LPA) (magnitude times 
2); + Time Difference Error (TDE); — Cumulative TDE. 



79 



5. KODIAK 




(a) Kodiak Station Zulu 



Figure 35. CALOC Strip Chart for Kodiak Loran-C Chain: o Linear Phase Adjust- 
ment (LPA) (magnitude times 2); -f- Time Difference Error (TDE); — Cumulative 
TDE. 
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APPENDIX B. STRIP CHARTS GENERATED 
BY PID CONTROLLER 

1. SOUTH EAST UNITED STATES (SEUS) 
a. 18 January 1997 




(a) SEUS Station Yankee 




(b) SEUS Station Zulu 

Figure 36. PID Controller Strip Chart for South East United States (SEUS) Loran-C 
Chain recorded on 18 January 1997: o Linear Phase Adjustment (LPA) (magnitude 
times 2); + Time Difference Error (TDE); — Cumulative TDE. 
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b. 20 January 1997 




(a) SEUS Station Whiskey 




(b) SEUS Station Xray 



82 



TDE (ns) TDE (ns) 




(c) SEUS Station Yankee 




(d) SEUS Station Zulu 

Figure 37. PID Controller Strip Chart for South East United States (SEUS) Loran-C 
Chain recorded on 20 January 1997: o Linear Phase Adjustment (LPA) (magnitude 
times 2); + Time Difference Error (TDE); — Cumulative TDE. 
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2. SOUTH CENTRAL UNITED STATES (SOCUS) 
a. 18 January 1997 




(a) SOCUS Station Victor 




(b) SOCUS Station Whiskey 
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(c) SOCUS Station Xray 




(d) SOCUS Station Yankee 




(e) SOCUS Station Zulu 



Figure 38. PID Controller Strip Chart for South Central United States (SOCUS) 

Loran-C Chain recorded on 18 January 1997: o Linear Phase Adjustment (LPA) 

(magnitude times 2); + Time Difference Error (TDE); — Cumulative TDE. 
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b. 20 January 1997 




(a) SOCUS Station Victor 




(b) SOCUS Station Whiskey 
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(c) SOCUS Station Xray 




(d) SOCUS Station Yankee 




(e) SOCUS Station Zulu 

Figure 39. PID Controller Strip Chart for South Central United States (SOCUS) 
Loran-C Chain recorded on 20 January 1997: o Linear Phase Adjustment (LPA) 
(magnitude times 2); + Time Difference Error (TDE); — Cumulative TDE. 
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3. NORTH CENTRAL UNITED STATES (NOCUS) 
a. 4 May 1997 




(a) NOCUS Station Whiskey 




(b) NOCUS Station Xray 
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(c) NOCUS Station Yankee 



Figure 40. PID Controller Strip Chart for North Central United States (NOCUS) 
Loran-C Chain recorded on 4 May 1997: o Linear Phase Adjustment (LPA) (magni- 
tude times 2); -j- Time Difference Error (TDE); — Cumulative TDE. 
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b. 5 May 1997 




(a) NOCUS Station Whiskey 




(b) NOCUS Station Xray 
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(c) NOCUS Station Yankee 



Figure 41. PID Controller Strip Chart for North Central United States (NOCUS) 
Loran-C Chain recorded on 5 May 1997: o Linear Phase Adjustment (LPA) (magni- 
tude times 2); + Time Difference Error (TDE); — Cumulative TDE. 
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4. UNITED STATES WEST COAST (USWC) 
a. 4 May 1997 




(a) USWC Station Whiskey 




(b) USWC Station Xray 
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(c) USWC Station Yankee 

Figure 42. PID Controller Strip Chart for United States West Coast (USWC) Loran- 
C Chain recorded on 4 May 1997: o Linear Phase Adjustment (LPA) (magnitude 
times 2); + Time Difference Error (TDE); — Cumulative TDE. 
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b. 5 May 1997 




(a) USWC Station Whiskey 




(b) USWC Station Xray 
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(c) USWC Station Yankee 

Figure 43. PID Controller Strip Chart for United States West Coast (USWC) Loran- 
C Chain recorded on 5 May 1997: o Linear Phase Adjustment (LPA) (magnitude 
times 2); + Time Difference Error (TDE); — Cumulative TDE. 
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5. KODIAK 




(a) Kodiak Station Zulu 

Figure 44. PID Controller Strip Chart for Kodiak Loran-C Chain; o Linear Phase 
Adjustment (LPA) (magnitude times 2); + Time Difference Error (TDE); — Cumu- 
lative TDE. 
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APPENDIX C. STRIP CHARTS GENERATED 

BY KALMAN FILTER 

1. SOUTH EAST UNITED STATES (SEUS) 
a. 18 January 1997 




(a) SEUS Station Yankee 




(b) SEUS Station Zulu 



Figure 45. Kalman Filter Strip Chart for South East United States (SEUS) Loran-C 
Chain recorded on 18 January 1997: o Linear Phase Adjustment (LPA) (magnitude 
times 2); + Time Difference Error (TDE); — Cumulative TDE. 
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b. 20 January 1997 




(a) SEUS Station Whiskey 




(b) SEUS Station Xray 
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(c) SEUS Station Yankee 




(d) SEUS Station Zulu 



Figure 46. Kalman Filter Strip Chart for South East United States (SEUS) Loran-C 
Chain recorded on 20 January 1997: o Linear Phase Adjustment (LPA) (magnitude 
times 2); + Time Difference Error (TDE); — Cumulative TDE. 
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2. SOUTH CENTRAL UNITED STATES (SOCUS) 
a. 18 January 1997 




(a) SOCUS Station Victor 




(b) SOCUS Station Whiskey 
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(c) SOCUS Station Xray 




(d) SOCUS Station Yankee 




(e) SOCUS Station Zulu 



Figure 47. Kalman Filter Strip Chart for South Central United States (SOCUS) 
Loran-C Chain recorded on 18 January 1997: o Linear Phase Adjustment (LPA) 
(magnitude times 2); J- Time Difference Error (TDE); — Cumulative TDE. 
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b. 20 January 1997 




(a) SOCUS Station Victor 




(b) SOCUS Station Whiskey 
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(e) SOCUS Station Zulu 



Figure 48. Kalman Filter Strip Chart for South Central United States (SOCUS) 
Loran-C Chain recorded on 20 January 1997: o Linear Phase Adjustment (LPA) 
(magnitude times 2); + Time Difference Error (TDE); — Cumulative TDE. 
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3. NORTH CENTRAL UNITED STATES (NOCUS) 
a. 4 May 1997 
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(a) NOCUS Station Whiskey 




(b) NOCUS Station Xray 
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(c) NOCUS Station Yankee 



Figure 49. Kalman Filter Strip Chart for North Central United States (NOCUS) 
Loran-C Chain recorded on 4 May 1997: o Linear Phase Adjustment (LPA) (magni- 
tude times 2); -1- Time Difference Error (TDE); — Cumulative TDE. 
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b. 5 May 1997 




(a) NOCUS Station Whiskey 




(b) NOCUS Station Xray 
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(c) NOCUS Station Yankee 



Figure 50. Kalman Filter Strip Chart for North Central United States (NOCUS) 
Loran-C Chain recorded on 5 May 1997: o Linear Phase Adjustment (LPA) (magni- 
tude times 2); + Time Difference Error (TDE); — Cumulative TDE. 
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4. UNITED STATES WEST COAST (USWC) 
a. 4 May 1997 




(a) USWC Station Whiskey 




(b) USWC Station Xray 
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(c) USWC Station Yankee 



Figure 51. Kalman Filter Strip Chart for United States West Coast (USWC) Loran-C 
Chain recorded on 4 May 1997: o Linear Phase Adjustment (LPA) (magnitude times 
2); + Time Difference Error (TDE); — Cumulative TDE. 
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(b) USWC Station Xray 
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(c) USWC Station Yankee 

Figure 52. Kalman Filter Strip Chart for United States West Coast (USWC) Loran-C 
Chain recorded on 5 May 1997; o Linear Phase Adjustment (LPA) (magnitude times 
2); + Time Difference Error (TDE); — Cumulative TDE. 
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5. KODIAK 




(a) Kodiak Station Zulu 



Figure 53. Kalman Filter Strip Chart for Kodiak Loran-C Chain: o Linear Phase 
Adjustment (LPA) (magnitude times 2); + Time Difference Error (TDE); — Cumu- 
lative TDE. 
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APPENDIX D. MATLAB CODE 



This appendix contains seven sections of Matlab code, and associated function 
and variable descriptions. The first section contains a listing of the data sets and 
major functions with a brief description of their purpose. Also contained in Section 1 
are the function variables along with their definition and functions in which they are 
used. Sections 3 and 4 contain the functions pid and kalman, respectively, as well 
as their associated driver functions. Driver functions are script files used to make 
repeated calls to a function using numerous different input options without having to 
duplicate the basic function. Section 5 contains the Matlab code to preprocess the 
recorded data, and Section 6 contains the code used to generate the CALOC strip 
charts based on the recorded data. Finally, the code used to generate the plot titles 
is contained in Section 7. 



1 . 



DATA SETS, FUNCTION AND VARIABLE 
DESCRIPTIONS 



a. Data Sets 



kodiak.dat North Pacific (NORPAC) Station Zulu recorded on 11 Septem- 
ber 1996. 

nocus4my.dat North Central United States (NOCUS) Stations Whiskey, Xray 
and Yankee recorded on 4 May 1997. 

nocus5my.dat North Central United States (NOCUS) Stations Whiskey, Xray 
and Yankee recorded on 5 May 1997. 



seusl8jn.dat South East United States (SEUS) Stations Yankee and Zulu 
recorded on 18 January 1997. 

seus20jn.dat South East United States (SEUS) Stations Whiskey, Xray, 
Yankee and Zulu recorded on 20 January 1997. 
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socusl8j.dat 


South Central United States (SOCUS) Stations Victor, 
Whiskey, Xray, Yankee and Zulu recorded on 18 January 1997. 


socus20j.dat 


South Central United States (SOCUS) Stations Victor, 
Whiskey, Xray, Yankee and Zulu recorded on 20 January 1997. 


uswc4my.dat 


South East United States (SEUS) Stations Whiskey, Xray and 
Yankee recorded on 4 May 1997. 


uswc5my.dat 


South East United States (SEUS) Stations Whiskey, Xray and 
Yankee recorded on 5 May 1997. 



b. Functions 



firham.m 


This function returns the difference equation coefficients for a 
linear phase FIR filter multiplied by a Hamming window of 
length M and digital cut-off frequency, o>c- 


kalman.m 


This function implements a Kalman Filter to control the Time 
Difference Error (TDE_data) of the LORAN-C System. The 
two controlled parameters are the cumulative TDE.data (long- 
term error) and the current TDE.data (short-term error). The 
estimate for the current TDE.data is generated using a Kalman 
estimator, which is then used as one of the feedback parameters 
for the Kalman filter. The cumulative TDE.data is computed 
by adding the number of minutes divided by 60 times the cur- 
rent TDE.data to account for the number of averages per hour. 



parse-title. m This function parses the data file name string and generates a 
string for the title of the output graph. 



114 



pid.m 


This function implements a PID controller to control the Time 
Difference Error (TDE) of the LORAN-C System. The two 
controlled parameters are the cumulative TDE (long-term er- 
ror) and the current TDE (short-term error). The cumulative 
TDE is computed by adding the number of minutes divided 
by 60 times the current TDE to account for the number of 
averages per hour. 



c. Data Vector Variables 

seuslSjnY Data vector for SEUS Yankee on 18 January 1997. 



seuslSjnZ 


Data vector for SEUS Zulu on 18 January 1997. 


seus20jnW 


Data vector for SEUS Victor on 20 January 1997. 


seus20jnX 


Data vector for SEUS Whiskey on 20 January 1997. 


seus20jnY 


Data vector for SEUS Yankee on 20 January 1997. 


seus20jnZ 


Data vector for SEUS Zulu on 20 January 1997. 


socuslSjV 


Data vector for SOCUS Victor on 18 January 1997. 


socuslSjW 


Data vector for SOCUS Whiskey on 18 January 1997. 


socuslSjX 


Data vector for SOCUS Xray on 18 January 1997. 


socuslSjY 


Data vector for SOCUS Yankee on 18 January 1997. 


socuslSjZ 


Data vector for SOCUS Zulu on 18 January 1997. 


socus20jV 


Data vector for SOCUS Victor on 20 January 1997. 


socus20jW 


Data vector for SOCUS Whiskey on 20 January 1997. 


socus20jX 


Data vector for SOCUS Xray on 20 January 1997. 


socus20jY 


Data vector for SOCUS Yankee on 20 January 1997. 


socus20jZ 


Data vector for SOCUS Zulu on 20 January 1997. 


uswc4myW 


Data vector for USWC Whiskey on 4 May 1997. 


uswc4myX 


Data vector for USWC Xray on 4 May 1997. 
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uswc4myY 

uswc5myW 

uswc5myX 

uswcSmyY 

nocus4myW 

nocus4myX 

nocus4myY 

nocusSmyW 

nocus5myX 

nocus5myY 

kodiakZ 



Data vector for 
Data vector for 
Data vector for 
Data vector for 

Data vector for 
Data vector for 
Data vector for 
Data vector for 
Data vector for 
Data vector for 

Data vector for 



USWC Yankee on 4 May 1997. 

USWC Whiskey on 5 May 1997. 
USWC Xray on 5 May 1997. 

USWC Yankee on 5 May 1997. 

NOCUS Whiskey on 4 May 1997 
NOCUS Xray on 4 May 1997 
NOCUS Yankee on 4 May 1997 
NOCUS Whiskey on 5 May 1997 
NOCUS Xray on 5 May 1997 
NOCUS Yankee on 5 May 1997 

NORPAC Zulu on 11 September 1996 



d. Other Variables 



A Weight in the Matlab function DLQR con- 

straint equation associated with the feed- 
back parameters. 

B Weight in the Matlab function DLQR con- 

straint equation associated with the LPA. 



conv.data Result of the convolution between the FIR 
filter and data vector. 

corr20 Constant value used for applying corrections 

to the data to remove CALOC influence due 
to a 20 ns LPA. 



kalman 



kalman 



kalman, pid 



Data Correction 
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Data Correction 



corr40 

cum_TDE 

data.cols 

dec-data 

delay 

firJilter 
graph -title 

Kp 

Ki 

Kd 

L 

Ipa 



Constant value used for applying corrections 
to the data to remove CALOC influence due 
to a 40 ns LPA. 

Cumulative time difference error. 

Number of data points in the vector after 
decimation by factor of timeJnt. 

Reshaped result of convolution in prepara- 
tion for decimation. 

Constant value of actual delay experienced 
between an LPA being ordered and its effect 
being seen at the monitor station. 

DiflPerence equation coefficients for the linear 
phase FIR filter with Hamming window. 

String used to generate title for plot of strip 
chart. 

Proportional weight for PID Controller. 

Integral weight for PID Controller. 

Derivative weight for PID Controller. 

Control vector used to weight feedback pa- 
rameters TDE_est and cum_TDE. 

Locations of LPAs used for plotting the 
CALOC Strip Charts. 



kalman, pid 



kalman, pid 



kalman, pid 



Data Correction 



kalman, pid 



kalman, pid 



pid 



pid 



pid 



kalman 



kalman, pid 
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LPA 


Linear phase adjustment that is quantizer 
output. 


kalman, pid 


LPAJnt 


Constant value of the linear phase adjust- 
ment increment. 


kalman, pid 


minutes 


Number of minutes between TDE averages 
used to determine decimation interval. 


kalman, pid 


num.cols 


Number of columns in data set. 


kalman, pid 


num_rows 


Number of rows in data set. 


kalman, pid 


out 


Output of PID Controller or Kalman Esti- 
mator as input to quantizer. 


kalman, pid 


P 


Correlation matrix of the state error equa- 
tion. 


kalman 


plotlpa 


Vector of NaN used to generate LPA vector 
by indexing on quantizer output. 


kalman, pid 


Q 


Covariance of the system noise term associ- 
ated with the LPA in the Kalman model. 


kalman 


QJqr 


Weight in the Matlab function DLQR cost 
function associated with the feedback pa- 
rameters. 


kalman 


quant _var 


Half-step size for LPA quantizer. 


kalman, pid 


R 


Covariance of the sensor noise term associ- 
ated with the TDE in the Kalman model. 


kalman 


RJqr 


Weight in the Matlab function DLQR cost 
function associated with the LPA. 


kalman 
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kalman, pid 



TDE-data 

TDE_est 

timeJnt 

xaxis 



Data decimated by factor of timeJnt. 

Time difference estimate for Kalman filter. 

Number of data points that make up a dec- 
imation interval. 

Vector used to plot output values to com- 
mon reference in time. 



kalman 



kalman, pid 



kalman, pid 
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2. LINEAR PHASE FIR FILTER WITH HAMMING 
WINDOW 



function h = firham(M, omega_c) 
y, function h = firham(M, oraega_c) 

y. 

y, FIRHAM Returns the difference equation coefficients for a linear 
y, phase FIR filter multiplied by a Hamming window of length 

y, M and digital cutoff frequency, omega_c. 

y. Initialize variables for filter length and coefficients 
f ilter.length = 

diffeq_coeff = zeros (size(f ilter_length) ) ; 

y. Generate difference equation coefficients of linear phase FIR 
y. filter using equation (8.1.20) on page 584 from Proakis 
y. Digital Signal Processing 

diffeq.coeff = sin(omega_c * (f ilter_length - (M-l)/2)) ./ ... 

(pi * (f ilter_length - (M-l)/2)); 

y. If filter length is odd, then h(n) at (M-D/2 is omega_c/pi 
if rem(M,2)==l 

diffeq_coeff ((M-D/2 + 1) = omega_c/pi; 

end 

y. Generate Hamming window of length M 
ham_M = hamming (M) ; 

y, Multiply linear phase FIR filter by Heimming window of length M 
fir_ham = diffeq.coeff .* ham_M’; 

y. Make filter unit gain 
h = f ir_ham. /sum(f ir_ham) ; 
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3. PID CONTROLLER FUNCTION AND DRIVERS 



a. Function 



function out = pid(data_set , init_cond, graph_title) 
y. LORAN-C Thesis 
% 

y, PID This function implements a PID controller to control the Time 
y, time Difference Error (TDE) of the LORAN-C System. The two 

y, controlled parameters are the cumulative TDE (long-term error) 

% and the current TDE (short-term error) . The cumulative TDE is 

y, computed by adding the number of minutes divided by 60 times 

y, the current TDE to account for the number of averages per hour. 

y, 

y The inputs to the function are an integer data vector, the initial 
y, condition in integer form, and a string that will be the title of 
y, of the graph output . The minimum input required is a valid data 
y, vector. The initial condition or graph title may be omitted. If 
% either areomitted, their default values will be used. The default 
y, values for the other inputs are: 

% init_cond = 0; 

y. graph.title = ’ ’ ; 

y, Function call error checking. 

y Verify that function call has 1, 2 or 3 inputs, 
if nargin > 3 

error( ’There are too many arguments in the function call.’) 
elseif nargin == 0 

error (’ Invalid function call. At a minimiim, function call must 
have a data vector.’) 
end 

y Verify first entry is a valid integer data vector, 
if isstr(data_set) 

error (’Data set has been entered as a string. Data set must be a 
vector. ’ ) 
end 
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y, Verify each input is of the correct form. 

% 

y. 

if nargin == 3 

if isstr(init_cond) 

error (’ Initial condition entered as a string - must be an 
integer . 0 
else 

mbint (init_cond) 

end 

if "isstr(graph_title) 

errorC’Graph title is not a string.') 

end 

elseif nargin == 2 

if (isstr (init_cond) ) 

graph_title = init_cond; 
init_cond = 0; 

disp(’ Assuming second input to function is graph’) 

disp(’ title. Setting initial condition to 0.’) 

pause(5) 
else 

graph_title = ’ ' ; 

end 

else 

init_cond = 0; 
graph.title = ’ ’ ; 

end 

y. Generate linear phase FIR filter with Hamming window 
y, FIRHAM Returns the difference equation coefficients for a 
y, linear phase FIR filter multiplied by a Hamming window 

y, of length M and digital cutoff frequency, omega_c. 

fir_filter = firham(32, pi/45); 

y. Determine the size of the data file. 

[num_rows, num_cols] = size(data_set) ; 

y Constants for decimation. 

y. 

y. Determine number of data points after decimation. The length 
y. of time equating to the decimation interval is determined by the 
y, variable, minutes. Since each data point represents 10 seconds in 
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y, real time, there are (minutes x 6) data points in the decimation 

% interval. In this particular case, there are 45 data points in 

y, each 7.5 minute interval. Using the Matlab command floor rounds 

y. down the result of the operation to get an integer value used to 

% reshape the data for decimation. 

minutes = 7.5; 

time_int = minutes * 6; 

data_cols = floor (num_rows/time_int) ; 

y. Convolve FIR filter with data, 
conv.data = conv(f ir_f ilter, data.set) ; 

% Reshape data for decimation. 

dec_data = reshape (conv_data(l : (time_int*data_cols)) , time_int, 

data_cols) ; 

y. Decimate by factor of time.int to obtain a vector of IDEs taken 
y, at a 7.5 minute interval. These values are the input to the control 
'/o algorithm. 

TDE_data = dec_data(time_int , :); 

y. Initialize variables 
LPA_int = 20; 
quant_var = 10 ; 

Kp = -.35; 

Ki = -.12; 

Kd = -.2; 

out = zeros (size (TDE_data) ) ; 

LPA = zeros (size (TDE_data) ) ; 
cum_TDE = zeros (size (TDE.dat a) ) ; 
cum_TDE(l) = init.cond; 

y, PID Controller 

y, 

% This for loop implements a PID controller. The loop starting value 
% of 2 is a constraint of Matlab requiring the index of a vector to 
*/o be a positive, non-zero integer. 

for time = 2;data_cols 

% Calculate the cumulative time difference error (cum_TDE) by 
y, adding the current TDE (TDE_data) to the previous value for 
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y, the cumulative TDE. The current TDE is scaled by a factor of 
y. minutes/60 to account for the number of corrections per hour. 
cum_TDE(time) = cum_TDE(time-l) + TDE_data(time) * (minutes/60) ; 

y, Output of the PID controller prior to queintization. The integral 
y, weight, Ki, is applied to the cumulative TDE; the proportional 

y, weight, Kp, is applied to the current TDE; and the derivative 

y, weight, Kd, is applied to the difference between the current TDE 

y. and the value for the current TDE from the previous time step. 

out(time) = (Ki * cum_TDE(time) ) + (Kp * TDE_data(time)) + ... 

(Kd * (TDE_data(time) - TDE_data(time - 1))); 



y, LPA Decision 

y. 

y. The output of the PID controller is quantized to interface 
y, with the Lorein-C transmitter timing controller which accepts 
y. corrections in 20 ns intervals. It is a mid-tread quantizer 
y, with equal sized steps of 20 ns. The step size may be varied 
y, by adjusting the quantization variable, qucint.var. The output 
y, of the quantizer is the linear phase adjustment (LPA) . 
if out(time) >= 17*quant_vcir 
LPA (time) = 180; 

elseif out (time) >= 15*quant_var 
LPA (time) = 160; 

elseif out (time) >= 13*quant_var 
LPA (time) = 140; 

elseif out (time) >= ll*quant_var 
LPA (time) = 120; 
elseif out (time) >= 9*quant_var 
LPA (time) = 100; 
elseif out (time) >= 7*quant_var 
LPA (time) = 80; 

elseif out (time) >= 5*quant_var 
LPA (time) = 60; 

elseif out (time) >= 3*quant_var 
LPA (time) = 40; 
elseif out (time) >= quant.var 
LPA (time) = 20; 

elseif out (time) >= -quant _var 
LPA (time) = 0; 

elseif out (time) >= -3*quant_var 
LPA (time) = -20; 
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elseif out (time) >= -5*quant_var 
LPA(time) = -40; 

elseif out (time) >= -7*quant_var 
LPA(time) = -60; 

elseif out (time) >= -9*quant_var 
LPA(time) = -80; 

elseif out (time) >= -ll*quant_var 
LPA(time) = -100; 
elseif out (time) >= -13*quant_var 
LPA(time) = -120; 
elseif out (time) >= -15*quant_var 
LPA(time) = -140; 
elseif out(time) >= -17*quaLnt_var 
LPA(time) = -160; 
elseif out (time) >= -19*quemt_var 
LPA(time) = -180; 

end 

7o Apply LPA to the remaining data set. This is a function of 
7o non-real time implementation. If implemented in real-time, 
y, the LPA would be applied to the transmitter and the shift in 
7, TDE would be seen in the data points. 

TDE_data(time :data_cols) = TDE_data(time;data_cols) + LPA(time); 



end 

7c Generate vector to plot LPA. Dividing by zero generates a 
7o vector of NaN so that only the LPA locations determined in 
7. the for loop below will be plotted, 
plotlpa = zeros(size(LPA))/0.0; 

7o This for loop determines the locations of the LPAs by indexing 
7c the capability of Matlab to index a vector based on its non-zero 
7c values. The location of each LPA is transferred to the NaN 
7c vector generated above. Since Matlab does not plot NaN, the only 
7c values that are plotted are the LPA locations . 
for index=l :data_cols 
if LPA ( index) "=0 

plotlpa(index)=LPA(index) ; 

end 

end 
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*/, Generate the LORAN Strip Charts 

7 . 

% The plots are generated using subplot so that they appear with 
7 . similar resolution to the existing strip charts. The plots 
7 . include the TDE averages (TDE.data) , the cumulative TDE, and the 
7 . locations of the LPAs input by CALOC. The LPAs are plotted 
% with a multiplicative factor of 2 solely for ease of display. 

% The variable xaxis is used for common reference when plotting, 
xaxis = [1 :num_rows/time_int] /(60/minutes) ; 

figure 

subplot (311) , plot(xaxis, TDE_data, ’w+’, xaxis, cum_TDE(l :data_cols) , 

'w-’, xaxis, plotlpa*2, ’wo') 
xlabeK’Time (hours)’), ylabeK’TDE (ns)’) 
title (graph_title) 

2Lxis([0 24 -80 80]), grid on 



b. Drivers 



datal8seuscorr 

graph_title = parse_title(’seus_18jnY’) ; 
pid(seusl8jnY , -31, graph.title) ; 

graph_title = parse_title( ’seus_18jnZ’) ; 
pid(seusl8jnZ, -77, graph_title) ; 



data20seuscorr 

graph_title = parse_title( ’seus_20jnW’ ) ; 
pid(seus20jnW, -40, graph_title) ; 

graph_title = parse_title( ’seus_20jnX’ ) ; 
pid(seus20jnX, -60, graph_title) ; 

graph_title = parse_title(’seus_20jnY’) ; 
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pid(seus20jnY, -50, graph.title) ; 

graph.title = parse_title(’seus_20jnZ’) ; 
pid(seus20jnZ, -92, graph.title) ; 



datalSsocuscorr 

graph_title = parse_title(’socusl8jnV’) ; 
pid(socusl8jV , -55, graph.title) ; 

graph.title = parse_title(’socusl8jnW’) ; 
pid(socusl8jW, -17, graph.title) ; 

graph_title = parse_title(’socusl8jnX’) ; 
pid(socusl8jX, 21, graph_title) ; 

graph.title = parse_title(’socusl8jnY’) ; 
pid(socusl8jY, 16, graph.title) ; 

graph_title = parse_title(’socusl8jnZ’) ; 
pid(socusl8jZ, -20, graph_title) ; 



data20socuscorr 

graph.title = parse_title('socus20jnVO ; 
pid(socus20jV, -42, graph.title) ; 

graph_title = parse_title(’socus20jnW’) ; 
pid(socus20jW, -20, graph_title) ; 

graph.title = parse_title(’socus20jnXO ; 
pid(socus20jX, 20, graph_title) ; 

graph.title = parse_title(’socus20jnY’) ; 
pid(socus20jY, -42, graph.title) ; 

graph_title = parse_title(’socus20jnZ’) ; 
pid(socus20jZ, 4, graph_title) ; 
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data4nocuscorr 



graph_title = parse_title( ’nocus04myW' ) ; 
pid(nocus4myW, -35, graph_title) ; 

graph.title = parse_title( ’nocus04myX’) ; 
pid(nocus4myX, 30, graph_title) ; 

graph_title = parse_title(’nocus04myY’) ; 
pid(nocus4myY, 25, graph_title) ; 



data5nocuscorr 

graph_title = parse_title( ’nocus05myW’) ; 
pid(nocus5myW, -15, graph.title) ; 

graph.title = parse_title( ’nocus05myX’) ; 
pid(nocus5myX , -25, graph_title) ; 

graph_title = parse_title ( ’nocus05myY’ ) ; 
pid(nocus5myY, 5, graph_title) ; 



data4uswccorr 

graph_title = parse_title(’uswc_04myW’) ; 
pid(uswc4myW, 80, graph_title) ; 

graph_title = parse_title('uswc_04myX') ; 
pid(uswc4myX, -30, graph_title) ; 

graph_title = parse_title('uswc_04myY’) ; 
pid(uswc4myY, -40, graph_title) ; 
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dataSuswccorr 



graph_title = parse_title( ’uswc_05myWO ; 
pid(uswc5myW, 60, graph_title) ; 



graph.title = parse_title(’uswc_05myX’) ; 
pid(uswc5myX, 25, graph.title) ; 



graph_title = parse_title('uswc_05myY’) ; 
pid(uswc5myY, 25, graph_title) ; 



datakodiaJccorr 

graph_title = parse_title( ’kodikllspZ’) ; 
pid(kodiakZ, -15, graph_title) ; 
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4. KALMAN FILTER FUNCTION AND DRIVERS 



a. Function 



function out = kalman(data_set , init_cond, graph.title) 
% LORAN-C Thesis 

y. 



y. KALMAN 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 



This function implements a Kalmain Filter to control the 
Time Difference Error (TDE_data) of the LQRAN-C System. 

The two controlled parameters cire the cumulative TDE_data 
(long-term error) and the current TDE_data (short-term 
error) . The estimate for the current TDE_data is generated 
using a Kalmain estimator, which is then used as one of the 
feedback parameters for the Kalmcui filter. The cumulative 
TDE_data is computed by adding the number of minutes divided 
by 60 times the current TDE_data to account for the number 
of averages per hour. 



y. Function call error checking. 



y. Verify that function call has 1, 2 or 3 inputs, 
if nargin > 3 

error (’There are too mciny arguments in the function call.’) 
elseif nargin == 0 

error (’ Invalid function call. At a minimum, function call must 
have a data vector . ’ ) 
end 



y. Verify first entry is a valid integer data vector, 
if isstr (data_set) 

error (’Data set has been entered as a string. Data set must be a 
vector. ’) 
end 

y, Verify each input is of the correct form. 

y. 

y. 

if nargin == 3 

if isstr (init_cond) 
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error (’Initial condition entered as a string - must be an 
integer. ’ ) 
else 

mbint (init_cond) 

end 

if "isstr (graph_title) 

errorC’Graph title is not a string.’) 

end 

elseif naxgin == 2 

if (isstr (init_cond)) 

graph_title = init_cond; 
init_cond = 0; 

disp(’ Assuming second input to function is graph’) 

disp(’ title. Setting initial condition to 0.’) 

pause (5) 
else 

graph_title = ’ ’ ; 

end 

else 

init_cond = 0; 
graph_title = ’ ’ ; 

end 

y. Generate linear phase FIR filter with Hamming window 
y. FIRHAM Returns the difference equation coefficients for a 
y, linear phase FIR filter multiplied by a Hamming window 

y. of length M eind digital cutoff frequency, omega.c. 

fir_filter = firham(32, pi/45); 

y, Determine the size of the data file. 

[num_rows, num_cols] = size(data_set) ; 

y. Constants for decimation. 

y. 

y. Determine niimber of data points after decimation. The length 
y, of time equating to the decimation interval is determined by the 
y, variable, minutes. Since each data point represents 10 seconds in 
y. real time, there are (minutes x 6) data points in the decimation 
y. interval. In this particular case, there are 45 data points in 
y. each 7.5 minute interval. Using the Matlab command floor rounds 
y. down the result of the operation to get an integer value used to 
y. reshape the data for decimation. 
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minutes = 7.5; 

time_int = minutes * 6; 

data_cols = floor (num_rows/time_int) ; 

% Convolve FIR filter with data. 
conv_data = conv(f ir_f ilter , data_set) ; 

*/o Reshape data for decimation. 

dec_data = reshape (conv_data(l : (time_int*data_cols)) , time_int, 

data_cols) ; 

*/, Decimate by factor of time_int to obtain a vector of TDE.datas taken 
*/« at a 7.5 minute interval. These values are the input to the control 
*/. algorithm. 

TDE.data = dec_data(time_int , :); 

y. Initialize variables 
LPA_int = 20; 
quant _var = 25; 

R = median(std(dec_data) ) ; 

Q = .5; 

p = zeros(size(TDE_data)) ; 
p(l) = 10; 

LPA = zeros (size (TDE_data) ) ; 
cum.TDE = zeros (size (TDE_data) ) ; 
cum_TDE(l) = init_cond; 

TDE_est = zeros (size (TDE_data) ) ; 

% Linear Quadratic Regulator 

*/o L = DLQR(A, B, Q, R) calculates the optimal feedback gain matrix 
*/. K such that the feedback law 

7o u[n] = -Kx[n] 

7o minimizes the cost function 

7o J = Sum {x’Qx + u’Ru} 

7. subject to the constraint equation: 

7. x[n+l] = Ax[n] + Bu[n] 

A = [1 0; 1 1] ; 

B = [1 ; 0] ; 

R_lqr = 1 ; 

Q_lqr = diag([5 5]); 

L = dlqr(A, B, Q_lqr, R_lqr) ; 



132 



*/, Kalman estimator 
for time = 2:data_cols 



*/. Apply LPA to the remaining data set. This is a function of 
*/, non-real time implementation. If implemented in real-time, 
y, the LPA would be applied to the transmitter and the shift in 
7. TDE would be seen in the data points . 

TDE_data(time;data_cols) = TDE_data(time:data_cols) + LPA(time-l); 

7, Estimator for the current time difference error (TDE_est) . 
TDE_est(time) = TDE.est (time-1) + LPA(time-l) + ... 

(p(time-l)/(R + p(time-l))) * (TDE_data(time-l) - ... 

TDE_est (time-1) ) ; 

p(time) = p(time-l) + Q - (p(time-l) '‘2/ (R + p(time-l))); 

7 Calculate the cumulative time difference error (cum_TDE) by 
7 adding the current TDE (TDE_data) to the previous value for 
7 the cumulative TDE. The current TDE is scaled by a factor of 
7 minutes/60 to account for the number of corrections per hour. 
cum_TDE(time) = cum_TDE(time-l) + (TDE_data(time) * (minutes/60)); 

7 Output of the Kalmam estimator prior to quanitization. The two 
7 parameters that are fed back are the estimate of the current TDE 
7 (TDE_est) and the actual cumulative TDE (cum_TDE) 
out(tirae) = -L * [TDE_est(time) ; cum_TDE(time)] ; 

7 LPA Decision 
7 

7 The output of the Kalman estimator is quantized to interface 
7 with the Loran-C treinsmitter timing controller which accepts 
7 corrections in 20 ns intervals. It is a mid-tread quamtizer 
7 with equal sized steps of 20 ns. The step size may be varied 
7 by adjusting the quantization variable, quant _var. The output 
7 of the quantizer is the linear phase adjustment (LPA) . 
if out (time) >= 17*quant_var 
LPA(time) = 180; 

elseif out (time) >= 15*quant_var 
LPA (time) = 160; 

elseif out (time) >= 13*quant_var 
LPA (time) = 140; 

elseif out(time) >= ll*quant_var 
LPA (time) = 120; 
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elseif out (time) >= 
LPA(time) = 100; 
elseif out (time) >= 
LPA(time) = 80; 
elseif out (time) >= 
LPA(time) = 60; 
elseif out(time) >= 
LPA(time) = 40; 
elseif out (time) >= 
LPA(time) = 20; 
elseif out (time) >= 
LPA(time) = 0; 
elseif out (time) >= 
LPA(time) = -20; 
elseif out(time) >= 
LPA(time) = -40; 
elseif out (time) >= 
LPA(time) = -60; 
elseif out(time) >= 
LPA(time) = -80; 
elseif out (time) >= 
LPA(time) = -100 
elseif out (time) >= 
LPA(time) = -120 
elseif out (time) >= 
LPA(time) = -140 
elseif out (time) >= 
LPA(time) = -160 
elseif out (time) >= 
LPA(time) = -180 

end 



9*qucint_var 

7*qucint_var 

5*quaiit_var 

3*quant_var 

quant_var 

-quant _var 

-3*quant_var 

-5*quant_var 

-7*quant_var 

-9*quant_var 

-ll*qucuit_var 

-13*quant_var 

-15*quant_var 

-17* quant _var 

-19*quant_var 



end 

’/, Generate vector to plot LPA. Dividing by zero generates a 
y. vector of NaN so that only the LPA locations determined in 
y, the for loop below will be plotted, 
plotlpa = zeros (size (LPA) )/0.0; 

y. This for loop determines the locations of the LPAs by indexing 
y. the capability of Matlab to index a vector based on its non-zero 
y. values. The location of each LPA is transferred to the NaN 
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’/. vector generated above. Since Matlab does not plot NaN, the only 
% values that are plotted are the LPA locations, 
for index=l :data_cols 
if LPA ( index) ~=0 

plotlpa(index)=LPA (index) ; 

end 

end 

*/, Generate the LORAN Strip Charts 

y, The plots are generated using subplot so that they appear with 
y, similar resolution to the existing strip charts. The plots 
y, include the TDE averages (TDE_data) , the cumulative TDE, and the 
y. locations of the LPAs input by CALOC. The factor of 0.125 on 
y. the cumulative error is due to the frequency of the averages 
y (i.e., 7.5 minutes is 1/8 of aji hour). The LPAs are plotted 
y. with a multiplicative factor of 3 solely for ease of display, 
xaxis = [1 :num_rows/time_int] /(60/minutes) ; 

figure 

subplot(311) , plot(xaxis, TDE_data, ’w+’, xaxis, cum_TDE, ’w-’, 

xaxis, plotlpa*2, ’wo’) 
xlabeK’Time (hours)’), ylabeK’TDE (ns)’) 
title (graph_title) 
axis([0 24 -80 80]), grid on 

b. Drivers 



datal8seuscorr 

graph_title = parse_title( ’ seus_18jnY’ ) ; 
kalman(seusl8jnY , -31, graph_title) ; 

graph_title = parse_title(’seus_18jnZ’) ; 
kalman(seusl8jnZ, -77, graph.title) ; 
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data20seuscorr 



graph.title = parse_title(’seus_20jnW’) ; 
kalman(seus20jnW, -40, graph_title) ; 

graph_title = parse_title('seus_20jnX') ; 
kalman(seus20jnX, -60, graph_title) ; 

graph_title = parse_title(’seus_20jnY’) ; 
kalman(seus20jnY, -50, graph_title) ; 

graph.title = parse_title( 'seus_20jnZ') ; 
kalman(seus20jnZ, -92, graph_title) ; 



datalSsocuscorr 

graph.title = parse_title(’socusl8jnV’) ; 
kalman(socusl8jV, -55, graph_title) ; 

graph_title = parse_title(’socusl8jnW') ; 
kalman(socusl8jW, -17, graph_title) ; 

graph.title = parse_title('socusl8jnX’) ; 
kalman(socusl8jX , 21, graph_title) ; 

graph_title = parse_title( ’ socusl8jnY' ) ; 
kalmaii(socusl8jY, 16, graph_title) ; 

graph_title = parse_title('socusl8jnZ') ; 
kalman(socusl8jZ, -20, graph.title) ; 



data20socuscorr 

graph.title = parse_title(’socus20jnV’) ; 
kalmcm(socus20jV, -42, graph_title) ; 

graph_title = parse_title( ’socus20jnW’ ) ; 
kalman(socus20jW, -20, graph.title) ; 
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graph_title = parse_title(’socus20jnX’) ; 
kalmam(socus20jX, 20, graph_title) ; 

graph_title = parse_title(’socus20jnY’) ; 
kalman(socus20jY, -42, graph.title) ; 

graph_title = parse_title(’socus20jnZ’) ; 
kalmein(socus20jZ, 4, graph_title) ; 



data4nocuscorr 

graph_title = parse_title(’nocus04myW’) ; 
kalman(nocus4myW, -35, graph.title) ; 

graph.title = parse_title( ’nocus04myX’) ; 
kalmein(nocus4rayX, 30, graph_title) ; 

graph_title = parse_title(’nocus04myY’) ; 
kalra£in(nocus4myY, 25, graph_title) ; 



data5nocuscorr 

graph.title = parse_title(’nocus05myW’) ; 
kalm2m(nocus5myW, -15, graph_title) ; 

graph.title = parse_title( ’nocus05rayX’) ; 
kalmeua(nocus5myX, -25, graph.title) ; 

graph_title = parse_title(’nocus05myY’) ; 
kalman(nocus5myY, 5, graph_title) ; 



data4uswccorr 

graph_title = parse_title( ’uswc_04myW’ ) ; 
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kalman(uswc4myW, 80, graph_title) ; 



graph_title = parse_title( ’uswc_04myX’) ; 
kalman(uswc4myX, -30, graph.title) ; 

graph_title = parse_title(’uswc_04myY’) ; 
kalman(uswc4myY, -40, graph_title) ; 



dataSuswccorr 

graph_title = parse_title( ’uswc_05myW’) ; 
kalman(uswc5myW, 60, graph_title) ; 

graph_title = parse_title(’uswc_05myXO ; 
kalman(uswc5myX, 25, graph_title) ; 

graph_title = parse_title( ’uswc_05myY’) ; 
kalmaii(uswc5myY, 25, graph_title) ; 



datakodiaJccorr 

graph_title = parse_title(’kodikllspZ’) ; 
kalman(kodiakZ , -15, graph_title) ; 
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5. DATA CORRECTION 



•/. LORAN-C Thesis 

7. 

7, datalSseuscorr .m 

7. 

7o Data Set: Malone-SEUS 18Jaii97 

7. Program removes influence of CALOC from data. The average time 
7. delay was determined by observing the absolute time delay of 
7. each CALOC-generated LPA and talking the mean. 2 sigma outliers 
7o were removed to obtain the final solution of 56 time steps, or 
7. 560 seconds (9+ minutes) . Each LPA was removed by applying the 
7o opposite sign and linearly time phasing in the correction to 
7o obtain uncorrected data. 

7o Load data set 
load seusl8jn.dat; 

7o Remove bias from data. Original data has a range of 0-99 which 
7o equates to the actual rainge of -50 to +49. Each single step 
7o increase in the data equates to an actual 10 ns increase in the 
7. TDE, therefore the data is scaled by 10. 
seusl8jn = (seusl8jn - 50) * 10; 

7o Load LORAN chain data vectors for Yankee and Zulu. 
seusl8jnY = seusl8jn( : ,8) ; 
seusl8jnZ = seusl8jn( : ,7) ; 

y. Determine the size of the data file. 

[num_rows, num_cols] = size(seusl8jn) ; 

7o Correction delay as determined by hand calculation of mean eind 
7o variance of average LPA time delays. The remainder of this 
7. section generates the vector of length 56 for a given LPA to 
7o linearly enter the correction, 
delay = 56; 

corr20 = (20/delay) .* [1 : delay] ’ ; 

7o Remove LPAs to get raw data. The locations of the corrections 
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20 ; 
+ 



’/. were determined by manually determining the beginning and 
"/, ending points of the individual LPA. The time-phased correction 
"/, was centered on the reinge of the actual delay for each LPA. 
seusl8jnY(667 :722) = seusl8jnY(667:722) + corr20; 
seusl8jnY(723:1120) = seusl8jnY(723; 1120) + 
seusl8jnY(1121:1176) = seusl8jnY(1121 : 1176) 
seusl8jnY(1177:2133) = seusl8jnY(1177;2133) 
seusl8jnY(2134:2189) = seusl8jnY(2134:2189) 
seusl8jnY(2190:3073) = seusl8jnY(2190:3073) 
seusl8jnY(3074:3129) = seusl8jnY (3074: 3129) 
seusl8jnY(3130:4018) = seusl8jnY(3130:4018) 
seusl8jnY(4019:4074) = seusl8jnY(4019:4074) 
seusl8jnY(4075:5373) = seusl8jnY(4075:5373) 
seusl8jnY (5374: 5429) = seusl8jnY(5374:5429) 
seusl8jnY(5430:5902) = seusl8jnY (5430: 5902) 
seusl8jnY(5903:5958) = seusl8jnY(5903:5958) 

= seusl8jnY(5959:6423) 

= seusl8jnY(6424:6479) 

= seusl8jnY(6480:7217) 

= seusl8jnY(7218:7273) 



seusl8jnY (5959: 6423) 
seusl8jnY(6424:6479) 
seusl8jnY(6480:7217) 
seusl8jnY(7218:7273) 



+ corr20; 
+ corr20; 



20 + corr20; 
40; 

40 
60 
60 
80; 

80 + corr20; 
100 ; 

100 - corr20; 
80 
80 
60 
60 
40 
40 



- corr20; 



- corr20; 



corr20; 

seusl8jnY(7274:num_rows) = seusl8jnY(7274 :num_rows) + 20; 



seusl8jnZ(2229 ; 
seusl8jnZ(2285 : 
seusl8jnZ(4005 : 
seusl8jnZ(4061 : 
seusl8jnZ(5371 ; 
seusl8jnZ(5427 : 
seusl8jnZ(5905 : 
seusl8jnZ(6345 : 
seusl8jnZ(6401 : 
seusl8jnZ(6690 
seusl8jnZ(6746 
seusl8jnZ(7425 
seusl8jnZ(7481 



2284) = seus 
4004) = seus 
4060) = seus 
5370) = seus 
5426) = seus 
5904) = seus 
5960) = seus 
6400) = seus 
6689) = seus 
6745) = seus 
7424) = seus 
7480) = seus 
num_rows) = 



18jnZ (2229: 2284) + 
18jnZ (2285: 4004) + 
18jnZ(4005:4060) + 
18jnZ(4061:5370) + 
18jnZ(5371:5426) + 
18jnZ(5427:5904) + 
18jnZ (5905: 5960) + 
18jnZ (6345: 6400) - 
18jnZ (6401: 6689) - 
18jnZ (6690: 6745) - 
18jnZ (6746: 7424) - 
18jnZ (7425: 7480) - 
seusl8jnZ(7481 :num. 



+ corr20; 
- corr20; 



corr20 ; 

20 ; 

20 
40 
40 
20 

20 - corr20; 
corr20 ; 

20 ; 

20 - corr20; 
40; 

40 - corr20; 
rows) - 60; 



7. LORAN-C Thesis 
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•/. 

7, data20seuscorr .m 

7. 

7. Data Set: Malone-SEUS 20Jan97 

7o Program removes influence of CALOC from data. The average time 
7o delay was determined by observing the absolute time delay of 
7o each CALOC-generated LPA and taking the mean. 2 sigma outliers 
7o were removed to obtain the final solution of 56 time steps, or 
7o 560 seconds (9+ minutes) . Each LPA was removed by applying the 
7. opposite sign and linearly time phasing in the correction to 
7. obtain uncorrected data. 

7. Load data set 
load seus20jn.dat; 

7. Remove bias from data. Original data has a range of 0-99 which 
7. equates to the actual range of -50 to +49. Each single step 
7o increase in the data equates to am actual 10 ns increase in the 
7o TDE, therefore the data is scaled by 10. 
seus20jn = (seus20jn - 50) * 10; 

7o Load LORAN chain data vectors for Whiskey, Xray, Yankee, and Zulu. 

seus20jnW = seus20jn( : , 1) ; 

seus20jnX = seus20jn( : ,2) ; 

seus20jnY = seus20jn( : ,8) ; 

seus20jnZ = seus20jn( : ,7) ; 

7. Determine the size of the data file. 

[num_rows, num_cols] = size(seus20jn) ; 

7» Correction delay as determined by hand calculation of mean and 
7 variance of average LPA time delays . The remainder of this 
7. section generates the vector of length 56 for a given LPA to 
7. linearly enter the correction, 
delay = 56; 

corr20 = (20/delay) .* [1 :delay] ’ ; 
corr40 = (40/delay) .* [1 : delay] ’ ; 

7o Remove LPAs to get raw data. The locations of the corrections 
7o were determined by manually determining the beginning and 
7 ending points of the individual LPA. The time-phased correction 
7 was centered on the range of the actual delay for each LPA. 



141 



seus20jnW(1294 : 1349) = seus20jnW( 1294: 1349) + corr20; 
seus20jnW(1350: 1952) = seus20jnW( 1350: 1952) + 20; 
seus20jnW(1953:2008) = seus20jnW( 1953: 2008) - corr20; 
seus20jnW(3961 :4016) = seus20jnW(3961 :4016) - corr20; 
seus20jnW(4017:4631) = seus20jnW(4017:4631) - 20; 
seus20jnW (4632: 4687) = seus20jnW(4632:4687) - 20 - corr20; 
seus20jnW(4688:5511) = seus20jnW(4688:5511) - 40; 
seus20jnW(5512 : 5567) = seus20jnW(5512 : 5567) - 40 - corr20; 
seus20jnW(5568:7172) - seus20jnW(5568:7172) - 60; 
seus20jnW (7173: 7228) = seus20jnW(7173:7228) - 60 - corr20; 
seus20jnW(7229 :num_rows) = seus20jnW(7229:num_rows) - 80; 



seus20jnX(980: 

seus20jnX(1036 

seus20jnX(1794 

seus20jnX(2311 

seus20jnX(2367 

seus20jnX(3053 

seus20jnX(4129 

seus20jnX(4185 

seus20jnX(4559 

seus20jnX(4791 

seus20jnX(4847 

seus20jnX(5441 

seus20jnX(5497 

seus20jnX(5609 

seus20jnX(5665 

seus20jnX(6349 

seus20jnX(6405 

seus20jnX(7346 

seus20jnX(7402 



1035) = seus20jnX(980: 1035) + corr20; 
:1793) = seus20jnX(1036:1793) 

:1849) = seus20jnX(1794:1849) 

-.2366) = seus20jnX(2311:2366) 

:3052) = seus20jnX(2367:3052) 

:3108) = seus20jnX(3053:3108) 

:4184) = seus20jnX(4129:4184) 

:4558) = seus20jnX(4185:4558) 

:4614) = seus20jnX(4559:4614) 

:4846) = seus20jnX(4791 :4846) 

:5440) = seus20jnX(4847:5440) 

:5496) = seus20jnX(5441:5496) 

:5608) = seus20jnX(5497:5608) 

:5664) = seus20jnX (5609: 5664) 

:6348) = seus20jnX(5665:6348) 

:6404) = seus20jnX (6349: 6404) 

:7345) = seus20jnX(6405:7345) 

:7401) = seus20jnX(7346:740l) 



20 ; 

corr20; 

corr20; 

20 ; 

corr20; 

corr20; 

20 ; 

corr20; 

corr20; 

20 ; 

20 - corr20; 
40; 

40 - corr20; 
60; 

60 - corr20; 
80; 

80 - corr20; 



:num_rows) = seus20jnX(7402:num_rows) - 100; 



seus20jnY(301 :356) = seus20jnY(301 :356) + corr20; 
seus20jnY(357:886) = seus20jnY(357:886) + 20; 
seus20jnY(887:942) = seus20jnY(887 :942) + 20 + corr20; 
seus20jnY(943 : 1547) = seus20j nY (943 : 1547) + 40; 
seus20jnY(1548: 1603) = seus20jnY( 1548: 1603) + 40 + corr20; 
seus20jnY(1604:3126) = seus20jnY(1604:3126) + 60; 
seus20jnY(3127:3182) = seus20jnY(3127:3182) + 60 + corr20; 
seus20jnY (3183 : 3880) = seus20jnY (3183: 3880) + 80; 
seus20jnY (3881 : 3936) = seus20jnY (3881: 3936) + 80 + corr20; 
seus20jnY (3937: 4305) = seus20jnY (3937: 4305) + 100; 
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seus20jnY (4306: 4361) = 
seus20jnY (4362: 4472) = 
seus20jnY (4473: 4528) = 
seus20jnY (4529: 4800) = 
seus20jnY (4801: 4856) = 
seus20jnY (4857: 4963) = 
seus20jnY(4964:5019) = 
seus20jnY (5020: 5444) = 
seus20jnY (5445: 5500) = 
seus20jnY (5532: 5587) = 
seus20jnY (5588: 6357) = 
seus20jnY (6358: 6413) = 
seus20jnY (6414; 7401) = 
seus20jnY (7402: 7457) = 
seus20jnY (7458:num_rows 



seus20jnY(4306:4361) + 
seus20jnY (4362: 4472) + 
seus20jnY (4473: 4528) + 
seus20jnY (4529: 4800) + 
seus20jnY (4801: 4856) + 
seus20jnY (4857: 4963) + 
seus20jnY(4964:5019) + 
seus20jnY(5020:5444) + 
seus20jnY (5445: 5500) + 
seus20jnY(5532:5587) - 
seus20jnY (5588: 6357) - 
seus20jnY (6358: 6413) - 
seus20jnY(6414 :7401) - 
seus20jnY (7402; 7457) - 
) = seus20jnY(7458 :niini. 



seus20jnZ(1229:1284) = 
seus20jnZ(1285:2167) = 
seus20jnZ(2168;2223) = 
seus20jnZ(2224:3881) = 
seus20jnZ (3882: 3937) = 
seus20jnZ (3938; 4399) = 
seus20jnZ (4400: 4455) = 
seus20jnZ (4647: 4702) = 
seus20jnZ (4703: 4897) = 
seus20jnZ (4898: 4953) = 
seus20jnZ (4954: 5460) = 
seus20jnZ(5461;5516) = 
seus20jnZ (5517: 5625) = 
seus20jnZ (5626: 5681) = 
seus20jnZ (5682: 5712) = 
seus20jnZ (5713: 5768) = 
seus20jnZ (5769 : 5886) = 
seus20jnZ (5887: 5942) = 
seus20jnZ (5943: 6267) = 
seus20jnZ (6268: 6323) = 
seus20jnZ (6324: 7332) = 
seus20jnZ (7333: 7388) = 
seus20 jnZ (7389 : num_rows 



seus20jnZ(1229 : 1284) + 
seus20jnZ(1285:2167) + 
seus20jnZ(2168:2223) + 
seus20jnZ(2224:3881) + 
seus20jnZ (3882; 3937) + 
seus20jnZ (3938 : 4399) + 
seus20jnZ (4400 : 4455) + 
seus20jnZ (4647: 4702) - 
seus20jnZ (4703: 4897) - 
seus20jnZ (4898: 4953) - 
seus20jnZ (4954; 5460) - 
seus20jnZ(5461 :5516) - 
seus20jnZ(5517;5625) - 
seus20jnZ(5626:5681) - 
seus20jnZ(5682:5712) - 
seus20jnZ(5713:5768) - 
seus20jnZ(5769 :5886) - 
seus20jnZ(5887;5942) - 
seus20jnZ(5943:6267) - 
seus20jnZ (6268: 6323) - 
seus20jnZ (6324: 7332) - 
seus20jnZ(7333:7388) - 
) = seus20jnZ(7389 :num. 



100 - corr20; 
80; 

80 - corr20; 
60; 

60 - corr20; 
40; 

40 - corr20; 
20 ; 

20 - corr20; 
corr20; 

20 ; 

20 - corr20; 
40; 

40 - corr20; 
rows) - 60; 

corr20; 

20 ; 

20 + corr20; 
40; 

40 - corr20; 
20 ; 

20 - corr20; 
corr20; 

20 ; 

20 - corr20; 
40; 

40 - corr20; 
60; 

60 - corr20; 
80; 

80 - corr20; 
100 ; 

100 - corr20; 
120; 

120 - corr20; 
140; 

140 - corr20; 
.rows) - 160; 



143 



y. LORAN-C Thesis 

y, 

y, datalSsocuscorr .m 

y. 

y, Data Set : Malone-SOCUS 18Jan97 

y, Program removes influence of CALOC from data. The average time 
y, delay was determined by observing the absolute time delay of 
y, each CALOC-generated LPA and taking the meam. 2 sigma outliers 
y, were removed to obtain the final solution of 56 time steps, or 
y, 560 seconds (9+ minutes) . Each LPA was removed by applying the 
y, opposite sign and linearly time phasing in the correction to 
y, obtain uncorrected data. 

y. Load data set 
load s ocus 1 8 j . dat ; 

y. Remove bias from data. Original data has a range of 0-99 which 
y, equates to the actual reinge of -50 to +49. Each single step 
y, increase in the data equates to an actual 10 ns increase in the 
y, TDE, therefore the data is scaled by 10. 
socusl8j = (socusl8j - 50) * 10; 

y. Load LORAN chain data vectors for Victor, Whiskey, Xray, Yankee, 
y, and Zulu. 

socusl8jV = socusl8j ( : , 1) ; 
socusl8jW = socusl8j ( : ,2) ; 
socusl8jX = socusl8j ( : ,4) ; 
socusl8jY = socusl8j ( : ,6) ; 
socusl8jZ = socusl8j ( : ,9) ; 

y. Determine the size of the data file. 

[num_rows, num_cols] = size(socusl8j) ; 

y. Correction delay as determined by hcind calculation of mean and 
y, variance of average LPA time delays. The remainder of this 
y, section generates the vector of length 56 for a given LPA to 
y, linearly enter the correction, 
delay = 56; 

corr20 = (20/delay) .* [1 : delay] ’ ; 
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corr40 = (40/delay) .* [1 .delay] ’ ; 

'/, Remove LPAs to get raw data. The locations of the corrections 
7. were determined by manually determining the beginning and 
7. ending points of the individual LPA. The time-phased correction 
7. was centered on the rainge of the actual delay for each LPA. 
socusl8jV(2667;2722) = socusl8jV(2667 :2722) - corr20; 

socusl8jV(2723;4179) = socusl8jV(2723;4179) - 20; 

socusl8jV(4180:4235) = socusl8jV(4180;4235) + corr20; 

socusl8jV(5923:5978) = socusl8jV(5923:5978) - corr20; 

socusl8jV(5979:6883) = socusl8jV (5979: 6883) - 20; 

socusl8jV(6884:6939) = socusl8jV (6884: 6939) - 20 - corr20; 

socusl8jV(6940:7799) = socus 18 jV (6940 : 7799) - 40; 

socusl8jV(7800:7855) = socusl8jV(7800:7855) - 40 - corr20; 

socusl8jV(7856:num_rows) = socus 18 jV (7856 :num_rows) - 60; 

socusl8jW(1736:1791) = socusl8jW(1736: 1791) + corr20; 

socusl8jW(1792:2850) = socusl8jW(1792:2850) + 20; 

socusl8jW(2851 :2906) = socusl8jW(2851:2906) - corr20; 

socusl8jW(5313:5368) = socusl8jW(5313: 5368) - corr20; 

socusl8jW(5369 : 5999) = socusl8jW(5369: 5999) - 20; 

socusl8jW(6000 : 6055) = socusl8jW(6000 : 6055) - 20 - corr20; 

socusl8jW(6056:6798) = socusl8jW(6056 : 6798) - 40; 

socusl8jW(6799 : 6854) = socusl8jW (6799 : 6854) - 40 - corr20; 

socusl8jW(6855:num_rows) = socusl8jW(6855:num_rows) - 60; 

socusl8jX(1046 : 1101) = socusl8jX(1046: 1101) - corr20; 
socusl8jX(1102 :num_rows) = socusl8jX(1102 :num_rows) - 20; 

socusl8jY(1121: 1176) = socusl8jY(1121 : 1176) - corr20; 
socusl8jY(1177:1317) = socusl8jY(1177: 1317) - 20; 
socusl8jY(1318:1373) = socusl8jY (1318: 1373) - 20 - corr20; 
socusl8jY(1374:1727) = socusl8jY (1374: 1727) - 40; 
socusl8jY ( 1728 : 1783) = socusl8jY ( 1728 : 1783) - 40 + corr20; 
socusl8jY(1784:5929) = socusl8jY(1784:5929) - 20; 
socusl8jY(5930:5985) = socusl8jY(5930:5985) - 20 + corr20; 

socusl8jZ(543 : 598) = socusl8jZ(543 : 598) - corr20; 
socusl8jZ(599 : 1821) = socusl8jZ(599: 1821) - 20; 
socusl8jZ(1825: 1880) = socusl8jZ(1825: 1880) - 20 - corr20; 
socusl8jZ(1881:2903) = socusl8jZ(1881 :2903) - 40; 
socusl8jZ(2904:2959) = socusl8jZ(2904:2959) - 40 - corr20; 
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socusl8jZ(2960:4298) = 
socusl8jZ(4299:4354) = 
socusl8jZ(4355:4985) = 
socusl8jZ(4986:5041) = 
socusl8jZ (5042: 5076) = 
socusl8jZ(5077:5132) = 
socusl8jZ(5133:5806) = 
socusl8jZ (5807 : 5862) = 
socusl8jZ(6418:6473) = 
socusl8jZ(6474:7311) = 
socusl8jZ(7312:7367) = 
socus 18 jZ (7368 : num_rows 



socusl8jZ(2960:4298) - 


60; 




socusl8jZ(4299:4354) - 


60 - 


corr20 


socusl8jZ(4355:4985) - 


80; 




socusl8jZ(4986:5041) - 


80 + 


corr20 


socusl8jZ(5042:5076) - 


60; 




socusl8jZ(5077:5132) - 


60 + 


O 

O 

to 

o 


socusl8jZ(5133:5806) - 


40; 




socusl8jZ(5807:5862) - 


40 + 


o 

o 

O 


socusl8jZ(6418:6473) - 


corr20; 


socusl8jZ(6474:7311) - 


20; 




socusl8jZ(7312:7367) - 


20 - 


corr20 


) = socusl8jZ(7368 ;num. 


.rows) 


- 40; 



7. LORAN-C Thesis 

7. 

7. data20socuscorr .ra 

7. 

7. Data Set : Malone-SOCUS 20 Jcin97 

7. Program removes influence of CALOC from data. The average time 
7. delay was determined by observing the absolute time delay of 
7. each CALOC-generated LPA and talking the mean. 2 sigma outliers 
7. were removed to obtain the final solution of 56 time steps, or 
7i 560 seconds (9+ minutes) . Each LPA was removed by applying the 
7. opposite sign and linearly time phasing in the correction to 
% obtain uncorrected data. 

7. Load data set 
load socus20j.dat; 

7o Remove bias from data. Original data has a range of 0-99 which 
7. equates to the actual reinge of -50 to +49 . Each single step 
7o increase in the data equates to an actual 10 ns increase in the 
7o TDE, therefore the data is scaled by 10. 
socus20j = (socus20j - 50) * 10; 

7o Load LORAN chain data vectors for Victor, Whiskey, Xray, Yankee, 
7o and Zulu. 
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socus20jV = socus20j ( ; , 1) ; 
socus20jW = socus20j ( : ,2) ; 
socus20jX = socus20j ( : ,4) ; 
socus20jY = socus20j ( : ,6) ; 
socus20jZ = socus20j ( : ,9) ; 

y. Determine the size of the data file. 

[num_rows, num_cols] = size(socus20j) ; 

y. Correction delay as determined by heind calculation of mecin amd 
y. variance of average LPA time delays. The remainder of this 
y, section generates the vector of length 56 for a given LPA to 
% linearly enter the correction, 
delay = 56; 

corr20 = (20/delay) .* [1 : delay] ' ; 
corr40 = (40/delay) .*[1: delay] ’ ; 

% Remove LPAs to get raw data. The locations of the corrections 
y, were determined by manually determining the beginning and 

ending points of the individual LPA. The time-phased correction 
% was centered on the range of the actual delay for each LPA. 
socus20jV(222:277) = socus20jV(222 : 277) - corr20; 
socus20jV(278:6029) = socus20jV(278 :6029) - 20; 
socus20jV (6030: 6085) = socus20jV(6030:6085) - 20 - corr20; 
socus20jV(6086:num_rows) = socus20jV(6086:num_rows) - 40; 

socus20jW(368 :423) = socus20jW(368:423) - corr20; 
socus20jW(424:617) = socus20jW(424 : 617) - 20; 
socus20jW(618:673) = socus20jW(618 : 673) - 20 - corr20; 
socus20jW(674:1375) = socus20jW (674: 1375) - 40; 
socus20jW(1376:1431) = socus20jW(1376: 1431) - 40 + corr20; 
socus20jW( 1432: 3061) = socus20jW(1432 : 3061) - 20; 
socus20jW(3062:3117) = socus20jW(3062 : 3117) - 20 + corr20; 
socus20jW(3916:3971) = socus20jW(3916 : 3971) - corr20; 
socus20jW(3972:6423) = socus20jW(3972:6423) - 20; 
socus20jW(6424:6479) = socus20jW (6424: 6479) - 20 - corr20; 
socus20jW(6480:num_rows) = socus20jW(6480:num_rows) - 40; 

socus20jX(2701:2756) = socus20jX(2701 :2756) - corr20; 
socus20jX(2757:num_rows) = socus20jX(2757:num_rows) - 20; 

socus20jY(2870 : 2925) = socus20jY(2870:2925) - corr20; 
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socus20jY(2926:5627) = socus20jY(2926:5627) - 20; 
socus20jY(5628 : 5683) = socus20jY(5628:5683) - 20 - corr20; 
socus20jY(5684:6762) = socus20jY(5684:6762) - 40; 
socus20jY(6763:6818) = socus20jY(6763:6818) - 40 + corr20; 
socus20jY(6819 :num_rows) = socus20jY(6819 :num_rows) - 20; 



socus20jZ(349 :404) = 
socus20jZ (405: 3458) = 
socus20jZ(3459:3514) 
socus20jZ(3515 :4083) 
socus20jZ(4084:4139) 
socus20jZ(4140:4691) 
socus20jZ (4692: 4747) 
socus20 j Z (4748 : 6450) 
socus20jZ(6451 :6506) 
socus20jZ(6507 :7089) 
socus20jZ(7090:7145) 



socus20jZ(349 :404) - corr20; 



socus20jZ(405 :3458) - 


20; 


socus20jZ (3459: 3514) 


- 20 - 


socus20jZ(3515:4083) 


- 40; 


socus20jZ(4084:4139) 


- 40 - 


socus20jZ(4140:4691) 


- 60; 


socus20jZ (4692: 4747) 


- 60 + 


socus20jZ(4748:6450) 


- 40; 


socus20jZ(6451 :6506) 


- 40 + 


socus20jZ(6507:7089) 


- 20; 


socus20jZ(7090:7145) 


- 20 + 



corr20; 
corr20 ; 
corr20 ; 
corr20; 
corr20; 



•/. LORAN-C Thesis 

•/. 

% data4nocuscorr .m 

•/. 

*/, Data Set : Middletown-NOCUS 4May97 

*/, Program removes influence of CALOC from data. The average time 
*/, delay was determined by observing the absolute time delay of 
y, each CALOC-generated LPA and taking the mean. 2 sigma outliers 
y were removed to obtain the final solution of 56 time steps, or 
y. 560 seconds (9+ minutes) . Each LPA was removed by applying the 
y. opposite sign and linearly time phasing in the correction to 
y, obtain uncorrected data. 

y. Load data set 
load nocus4my.dat; 

y. Remove bias from data. Original data has a range of 0-99 which 
y, equates to the actual range of -50 to +49. Each single step 
y, increase in the data equates to an actual 10 ns increase in the 
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*/. TDE, therefore the data is scaled by 10. 
nocus4my = (nocus4my - 50) * 10; 

y. Load LORAN chain data vectors for Whiskey, Xray, Yankee. 
nocus4myW = nocus4my ( : , 10) ; 
nocus4myX = nocus4my ( : , 12) ; 
nocus4myY = nocus4my ( : ,3) ; 

*/o Determine the size of the data file. 

[num_rows, num_cols] = size(nocus4my) ; 



y. Correction delay as determined by hand calculation of mean and 
y« varicince of average LPA time delays . The remainder of this 
y, section generates the vector of length 56 for a given LPA to 
y, linearly enter the correction, 
delay = 56; 

corr20 = (20/delay) .* [1 : delay] ' ; 
corr40 = (40/delay) .♦ [1 : delay] ’ ; 



y. Remove LPAs to get raw data. The locations of the corrections 
y, were determined by manually determining the beginning and 
y, ending points of the individual LPA. The time-phased correction 
y, was centered on the reinge of the actual delay for each LPA. 
nocus4myW(565 :620) = nocus4myW(565 : 620) - corr40; 
nocus4myW(621 : 1033) = nocus4myW (621 : 1033) - 40; 



nocus4myW(1034 : 1089) 
nocus4myW ( 1090 : 3507 ) 
nocus4myW (3508 : 3563) 
nocus4myW (3564 : 4844) 
nocus4myW (4845 : 4900) 
nocus4myW(4901 :5922) 
nocus4myW(5923 ; 5978) 
nocus4myW (5979 : 7039) 
nocus4myW(7040 : 7095) 
nocus4myW(7096 ;7951) 
nocus4myW (7952 : 8007 ) 



nocus4myW(1034; 1089) 
nocus4myW ( 1090 ; 3507 ) 
nocus4myW (3508 : 3563) 
nocus4myW (3564 : 4844) 
nocus4myW (4845 : 4900) 
nocus4myW (4901 : 5922) 
nocus4myW (5923 : 5978) 
nocus4myW (5979 : 7039) 
nocus4myW(7040 ; 7095) 
nocus4myW (7096 : 7951) 
nocus4myW (7952 ; 8007 ) 



40 - 
60; 
60 - 
80; 
80 + 
60; 
60 + 
40; 
40 + 
20 ; 
20 - 



corr20; 



corr20 ; 



corr20; 



corr20 ; 



corr20; 



corr20 ; 



nocus4myW(8008:num_rows) = nocus4myW(8008:num_rows) - 40; 



nocus4myX(2136:2191) = nocus4myX (2136 : 2191) - corr20; 
nocus4myX (2192: 6456) = nocus4myX (2192: 6456) - 20; 
nocus4myX (6457 : 6512) = nocus4myX (6457: 6512) - 20 + corr20; 
nocus4myX(7047 :7102) = nocus4myX (7047: 7102) + corr20; 
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nocus4myX(7103:7976) 

nocus4myX(7977:8032) 



nocus4myX(7103:7976) + 20; 
nocus4myX (7977 : 8032) + 20 - corr20; 



nocus4myY(1715;1770) 
nocus4myY(1771 : 2213) 
nocus4myY(2214:2269) 
nocus4myY (2270 : 2650) 
nocus4myY(2651 :2706) 
nocus4myY (2707: 5650) 
nocus4myY (5651 : 5706) 
nocus4myY (5707: 6525) 
nocus4myY(6526:6581) 
nocus4myY (6582 : 7035) 
nocus4myY (7036 : 7091) 
nocus4myY (7092: 7966) 
nocus4myY (7967: 8022) 



= nocus4myY(1715:1770) 
= nocus4myY(1771:2213) 
= nocus4myY(2214:2269) 
= nocus4myY(2270:2650) 
= nocus4myY(2651:2706) 
= nocus4myY(2707:5650) 
= nocus4myY(5651:5706) 
= nocus4myY(5707:6525) 
= nocus4myY(6526:6581) 
= nocus4myY (6582 -.7035) 
= nocus4myY (7036:7091) 
= nocus4myY(7092:7966) 
= nocus4myY (7967: 8022) 



corr20; 

20 ; 

20 - corr20; 
40; 

40 - corr40; 
80; 

80 + corr20; 
60; 

60 + corr20; 
40; 

40 + corr20; 
20 ; 

20 + corr20; 



7. LORAN-C Thesis 

7 . 

7 . data5nocuscorr .m 
7 . 

7 . Data Set: Middletown-NOCUS 5May97 

7 . Program removes influence of CALOC from data. The average time 
7. delay was determined by observing the absolute time delay of 
7 . each CALOC-generated LPA and tahing the mean. 2 sigma outliers 
7 . were removed to obtain the final solution of 56 time steps, or 
7 . 560 seconds (9+ minutes) . Each LPA was removed by applying the 
7. opposite sign and linearly time phasing in the correction to 
7 obtain uncorrected data. 

7 Load data set 
load nocus5my.dat; 

7 Remove bias from data. Original data has a range of 0-99 which 
7 equates to the actual range of -50 to +49. Each single step 
7 increase in the data equates to an actual 10 ns increase in the 
7 TDE, therefore the data is scaled by 10. 
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nocusSmy = (nocus5my - 50) * 10; 

7, Load LORAN chain data vectors for Whiskey, Xray, Yankee. 
nocusBmyW = nocus5my( : , 10) ; 
nocusSmyX = nocus5my( : ,12) ; 
nocusSmyY = nocus5my( : ,3) ; 

7. Determine the size of the data file. 

[num_rows, num_cols] = size(nocus5my) ; 



7 » Correction delay as determined by hand calculation of mectn and 
% variance of average LPA time delays. The remainder of this 
7 . section generates the vector of length 56 for a given LPA to 
7. linearly enter the correction, 
delay = 56; 

corr20 = (20/delay) .♦ [1 : delay] ’ ; 
corr40 = (40/delay) .♦ [1 : delay] ' ; 

7 « Remove LPAs to get raw data. The locations of the corrections 
7. were determined by manually determining the beginning and 
7 . ending points of the individual LPA. The time-phased correction 
7 . was centered on the range of the actual delay for each LPA. 
nocus5myW(6554:6609) = nocus5myW(6554;6609) + corr20; 
nocus5myW(6610:num_rows) = nocus5myW(6610 :num_rows) + 20; 

nocus5myX(523:578) = nocus5myX(523 ; 578) - corr20; 
nocus5myX (579: 4687) = nocus5myX (579 : 4687) - 20; 
nocus5myX (4688: 4743) = nocus5myX (4688: 4743) -20 + corr20; 



nocus5myY(530:585) = nocus5myY(530:585) - corr20; 
nocus5myY (586: 1309) = nocus5myY(586 : 1309) - 20; 



nocus5myY(1310: 
nocus5myY ( 1366 : 
nocus5myY (2214 : 
nocus5myY (2270 : 
nocus5myY (3302 : 
nocus5myY (3358 : 
nocus5myY(5125: 
nocus5myY(5181 : 



1365) = nocus5myY(1310:1365) 
2213) = nocus5myY(1366:2213) 
2269) = nocus5myY (2214: 2269) 
3301) = nocus5myY(2270:3301) 
3357) = nocus5myY (3302: 3357) 
5124) = nocus5myY (3358: 5124) 
5180) = nocus5myY(5125:5180) 
6560) = nocus5myY(5181:6560) 
6616) = nocus5myY (6561: 6616) 



20 - corr20; 
40; 

40 - corr20: 
60; 

60 - corr20: 
80; 

80 + 

60; 

60 + 



corr20; 



corr40; 



nocus5myY (6561 : 
nocus5myY(6617 :num_rows) = nocus5myY(6617 :num_rows) - 20; 
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•/. LORAN-C Thesis 

y. 

*/. data4uswccorr .m 

•/. 

*/, Data Set: Middiet own-USWC 4May97 

*/. Program removes influence of CALOC from data. The average time 
'/, delay was determined by observing the absolute time delay of 
7, each CALOC-generated LPA and taking the mean. 2 sigma outliers 
7. were removed to obtain the final solution of 56 time steps, or 
7o 560 seconds (9+ minutes) . Each LPA was removed by applying the 
7. opposite sign and linearly time phasing in the correction to 
7. obtain uncorrected data. 

7. Load data set 
load uswc4my.dat; 

7« Remove bias from data. Original data has a range of 0-99 which 
7o equates to the actual range of -50 to +49. Each single step 
7o increase in the data equates to an actual 10 ns increase in the 
7o TDE, therefore the data is scaled by 10. 
uswc4my = (uswc4my - 50) * 10; 

7. Load LORAN chain data vectors for Whiskey, Xray, Yankee. 
uswc4myW = uswc4my( : , 1) ; 
uswc4myX = uswc4my ( : , 10) ; 
uswc4myY = uswc4my(: , 12) ; 

7o Determine the size of the data file. 

[num_rows, num_cols] = size(uswc4my) ; 

7. Correction delay as determined by hand calculation of mean and 
7o variance of average LPA time delays. The remainder of this 
7. section generates the vector of length 56 for a given LPA to 
7o linearly enter the correction, 
delay = 56; 

corr20 = (20/delay) .* [1 ; delay] ’ ; 

7. Remove LPAs to get raw data. The locations of the corrections 
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"/. were determined by meinually determining the beginning eind 

y. ending points of the individual LPA. The time-phased correction 

y. was centered on the reuige of the actual delay for each LPA. 

uswc4myW(175:230) = uswc4myW(175 : 230) + corr20; 

uswc4myW(231:2027) = uswc4myW(231 : 2027) + 20; 

uswc4myW ( 2028 : 2083) = uswc4myW(2028:2083) + 20 + corr20; 

uswc4myW(2084;2205) = uswc4myW(2084:2205) + 40; 

uswc4myW(2206 : 2261) = uswc4myW(2206:2261) + 40 + corr20; 

uswc4myW(2262:3880) = uswc4myW(2262:3880) + 60; 

uswc4myW(3881 : 3936) = uswc4myW (3881 : 3936) + 60 + corr20; 

uswc4myW(3937:4203) = uswc4myW (3937 : 4203) + 80; 

uswc4myW (4204 : 4259) = uswc4myW (4204: 4259) + 80 - corr20; 

uswc4myW (4260: 4661) = uswc4myW (4260 : 4661) + 60; 

uswc4myW(4662:4717) = uswc4myW(4662:4717) + 60 - corr20; 

uswc4myW(4718:num_rows) = uswc4myW(4718 :num_rows) + 40; 

uswc4myX(1515 : 1570) = uswc4myX ( 1515 : 1570) - corr20; 
uswc4myX( 1571: 2200) = uswc4myX ( 1571 : 2200) - 20; 
uswc4myX (2201: 2256) = uswc4myX (2201 : 2256) - 20 - corr20; 
uswc4myX(2257:3076) = uswc4myX (2257 : 3076) - 40; 
uswc4myX (3077 : 3132) = uswc4myX (3077 : 3132) - 40 + corr20; 
uswc4myX(3133 : 6897) = uswc4myX (3133 : 6897) - 20; 
uswc4myX (6898: 6953) = uswc4myX (6898 : 6953) - 20 - corr20; 
uswc4myX(6954 :num_rows) = uswc4myX(6954:num_rows) - 40; 

uswc4myY(1751 : 1806) = uswc4myY ( 1751 : 1806) - corr20; 
uswc4myY( 1807: 2437) = uswc4myY( 1807: 2437) - 20; 
uswc4myY(2438:2493) = uswc4myY (2438: 2493) - 20 - corr20; 
uswc4myY (2494 -.3733) = uswc4myY (2494: 3733) - 40; 
uswc4myY(3734:3789) = uswc4myY (3734 : 3789) - 40 + corr20; 
uswc4myY(3790 :num_rows) = uswc4myY(3790:num_rows) - 20; 



y, LORAN-C Thesis 

y. 

y, data5uswccorr . m 

y. 

y. Data Set: Middletown-USWC 5May97 
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"/. Program removes influence of CALOC from data. The average time 
*/. delay was determined by observing the absolute time delay of 
*/. each CALOC-generated LPA and taking the mean. 2 sigma outliers 
’/, were removed to obtain the final solution of 56 time steps, or 
*/, 560 seconds (9+ minutes) . Each LPA was removed by applying the 
'/, opposite sign and linearly time phasing in the correction to 
*/. obtain uncorrected data. 



*/. Load data set 
load uswc5my.dat; 

’/, Remove bias from data. Original data has a range of 0-99 which 
'/. equates to the actual reinge of -50 to +49. Each single step 
*/, increase in the data equates to am actual 10 ns increase in the 
’/, TDE, therefore the data is scaled by 10. 
uswc5my = (uswc5my - 50) * 10; 

y. Load LORAN chain data vectors for Whiskey, Xray, Yankee. 
uswc5myW = uswc5my( : , 1) ; 
uswc5myX = uswc5my( : , 10) ; 
uswc5myY = uswc5my( : , 12) ; 

y. Determine the size of the data file. 

[num_rows, num_cols] = size(uswc5my) ; 



y, Correction delay as determined by hand calculation of mean and 
y, variance of average LPA time delays. The remainder of this 
% section generates the vector of length 56 for a given LPA to 
y. linearly enter the correction, 
delay = 56; 

corr20 = (20/delay) .* [1 : delay] ’ ; 
corr40 = (40/delay) .* [1 ; delay] ’ ; 



y. Remove LPAs to get raw data. The locations of the corrections 
y, were determined by manually determining the beginning and 
y, ending points of the individual LPA. The time-phased correction 
y. was centered on the range of the actual delay for each LPA. 
uswc5myW(2016:2071) = uswc5myW (2016 : 2071) + corr20; 

= uswc5myW (2072: 2726) 

= uswc5myW (2727: 2782) 

= uswc5myW (2783: 3366) 

= uswc5myW (3367 : 3422) 



uswc5myW(2072:2726) 
uswc5myW (2727: 2782) 
us wc5my W (2783:3366) 
uswc5myW (3367 : 3422) 



20 ; 

20 + corr20; 
40; 

40 + corr20; 
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uswc5myW(3423:4202) = uswc5myW(3423 ; 4202) + 60; 
uswcSmyW (4203 : 4258) = uswc5myW(4203:4258) + 60 - corr20; 
uswc5myW(4259:5529) = uswc5myW(4259 : 5529) +40; 
uswc5myW ( 5530 : 5585) = uswc5myW(5530:5585) +40 - corr20; 
uswc5myW(5586 :num_rows) = uswc5myW(5586:mim_rows) + 20; 

uswc5myX(1470 : 1525) = uswc5myX( 1470: 1525) - corr20; 
uswc5myX(1526:3341) = uswc5myX ( 1526 : 3341) - 20; 
uswc5myX (3342: 3397) = uswc5myX (3342 : 3397) - 20 + corr40; 
uswc5myX(3398:3610) = uswc5myX (3398: 3610) + 20; 
uswc5myX (361 1 : 3666) = uswc5myX (361 1 : 3666) + 20 - corr20; 
uswc5myX (4197: 4252) = uswc5myX (4197 : 4252.) - corr20; 
uswc5myX(4253 :num_rows) = uswc5myX(4253:num_rows) - 20; 

uswc5myY (2122: 2177) = uswc5myY(2122 : 2177) - corr20; 
uswc5myY (2178: 3345) = uswc5myY(2178:3345) - 20; 
uswc5myY (3346 : 3401) = uswc5myY(3346:3401) - 20 + corr40; 
uswc5myY (3402: 4174) = uswc5myY (3402 : 4174) + 20; 
uswc5myY(4175 : 4230) = uswc5myY(4175:4230) + 20 - corr20; 
uswc5myY (7021: 7076) = uswc5myY (7021 : 7076) - corr20; 
uswc5myY(7077:n\im_rows) = uswc5myY(7077:num_rows) - 20; 



•/. LORAN-C Thesis 

y. 

7, datakodiakcorr .m 

7 . 

7, Data Set : Kodiak Zulu 

7o Program removes influence of CALOC from data. The average time 
7, delay was determined by observing the absolute time delay of 
7, each CALOC-generated LPA eind taking the mean. 2 sigma outliers 
7, were removed to obtain the final solution of 56 time steps, or 
7, 560 seconds (9+ minutes) . Each LPA was removed by applying the 
7. opposite sign eind lineeirly time phasing in the correction to 
7. obtain uncorrected data. 

7o Load data set 
load kodiak.dat; 
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7, Remove bias from data. Original data has a rainge of 0-99 which 
7, equates to the actual range of -50 to +49. Each single step 
7o increase in the data equates to an actual 10 ns increase in the 
7 . TDE, therefore the data is scaled by 10. 
kodiak = (kodiak - 50) + 10; 

7 . Load LORAN chain data vector for Zulu. 
kodiakZ = kodiak( : , !)•; 



7 , Determine the size of the data file. 
[num_rows, num.cols] = size (kodiak) ; 



7 « Correction delay as determined by hand calculation of mean cuid 
7 « variance of average LPA time delays. The remainder of this 
7 « section generates the vector of length 56 for a given LPA to 
7o linearly enter the correction, 
delay = 56; 

corr20 = (20/delay) .* [1 .delay] ’ ; 



7 . Remove LPAs to get raw data. The locations of the corrections 
% were determined by manually determining the beginning ctnd 
7 . ending points of the individual LPA. The time-phased correction 
7 » was centered on the range of the actual delay for each LPA. 
kodiakZ (343: 398) = kodiakZ (343: 398) - corr20; 
kodiakZ(399: 1687) = kodiakZ (399 : 1687) - 20; 
kodiakZ (1688: 1743) = kodiakZ (1688: 1743) 

= kodiakZ(1981:2036) 

= kodiakZ (2037: 3619) 

= kodiakZ (3620: 3675) 



kodiakZ(1981:2036) 
kodiakZ(2037:3619) 
kodiakZ (3620: 3675) 



- 20 + corr20; 

- corr20; 

- 20 ; 

- 20 + corr20; 
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6. CALOC STRIP CHARTS 



•/. LORAN-C Thesis 
y. seuslS.m 

y. 

y, Data Set : Malone-SEUS 18Jan97 

y, Progreun generates LORAN Time Difference Error (TDE) vs Time strip 
y, charts similar to the output of the current CALOC strip charts, 
y. These strip charts will be used to compare the output of the new 
y, control algorithm to determine the increase in control accuracy. 

y. Clear stored variables 
clear 

y, Load data set 
load seusl8jn.dat; 

y. Remove bias from data. Original data has a range of 0-99 which 
y. equates to the actual range of -50 to +49. Each single step 
y. increase in the data equates to cin actual 10 ns increase in the 
y. TDE, therefore the data is scaled by 10. 
seusl8jn = (seusl8jn - 50) + 10; 

y. Load LORAN chain data vectors for Yankee and Zulu. 
seusl8jnY = seusl8jn( : ,8) ; 
seusl8jnZ = seusl8jn( : ,7) ; 

y. Determine the size of the data file. 

[num_rows, num_cols] = size (seusl8jn) ; 

*/o Determine number of data points after taking 7.5 minute averages, 
y. Since each data point represents 10 seconds in real time, there 
y, are 45 data points in each 7.5 minute average. Using the Matlab 
y. floor command rounds down the result of the operation to get an 
y, integer value to use in reshaping the data for averaging. The 
y. second line generates the vector for plotting the data. 
data_cols = floor (num_rowsy45) ; 
xaxis = [1 :num_rows/45] /8; 
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’/, Take 7.5 minute averages to generate data points for plotting 
7, similar to current strip charts. Averages are taken by reshaping 
7o the data into a matrix of dimensions data_cols x 45 and then tciking 
7, the mean in a column-wise fashion. This generates a vector of 
7. length data.cols of the 7.5 minute averages. 

avgY = meein(reshape(seusl8jnY(l : (45*data_cols)) , 45, data_cols)); 
avgZ = mecLn(reshape(seusl8jnZ(l ; (45*data_cols) ) , 45, data.cols)); 



7o After plotting the raw data, the locations of the LPAs were 
7o determined by comparing the generated plot with the actual strip 
7. charts. The locations are placed below to generate LPA vectors. 
Ipay = zeros(size(l :data_cols) ) ; 

Ipaz = zeros(size(l :data_cols) ) ; 



Ipay (16) = -20; 
Ipay (69) = -20; 
Ipay (120) = 20; 
Ipay (144) = 20; 



Ipay (26) = -20; 
Ipay (90) = -20; 
Ipay (132) = 20; 
Ipay (161) = 20; 



Ipay (48) = -20; 



Ipaz (50) = -20; 
Ipaz (120) = 20; 
Ipaz (149) = 20; 



Ipaz (90) = -20; 
lpaz(132) = 20; 
Ipaz (166) = 20; 



Ipaz (142) = 20; 



7o Generate vector to plot LPA. Dividing by zero generates a 
7o vector of NaN so that only the LPA locations determined in 
7o the for loop below will be plotted, 
plotlpay = zeros (size (Ipay) )/0 .0 ; 
plotlpaz = zeros(sizedpaz) )/0 .0 ; 



7. This for loop determines the locations of the LPAs by indexing 
7c the capability of Matlab to index a vector based on its non-zero 
7c values. The location of each LPA is treinsf erred to the NaN 
7c vector generated above. Since Matlab does not plot NaN, the only 
7c values that are plotted are the LPA locations, 
for index=l :data_cols 
if Ipay ( index) ~=0 

plotlpay (index) =lpay( index) ; 

end 

if Ipaz (index) ~=0 

plotlpaz ( index) =lpaz( index) ; 

end 
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end 



y, Generate the LORAN Strip Charts 

y. 

y, The plots cire generated using subplot so that they appear with 
y. similar resolution to the existing strip charts. The plots 
y, include the TDK averages (avgY) , the cumulative TDE, and the 
% locations of the LPAs input by CALOC. The initial value for 
y. the cumulative error is taken from the actual strip charts from 
y. CALOC. The factor of 0.125 on the cumulative error is due to 
% the frequency of the averages (i.e., 7.5 minutes is 1/8 of cin 
y, hour) . The LPAs are plotted with a multiplicative factor of 3 
y, solely for ease of display, 
f igure 

subplotOll) , plot(xaxis, avgY, ’w+') 
xlabelC’Time (hours)’), ylabel(’TDE (ns)’) 
hold on 

plot(xaxis, -31 + cumsum(avgY * 0.125), ’w-’) 
plot(xaxis, plotlpay * 3, ’wo’) 
grid on 

axis ( [0 24 -80 80] ) 

title(’SEUS Yainkee 18Jan - Plot of TDE vs Time (Strip Chart)’) 
figure 

subplot(311) , plot(x 2 Lxis, avgZ, ’w+’) 
xlabeK’Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, -77 + cumsum(avgZ * 0.125), ’w-’) 
pl6t(xaxis, plotlpaz * 3, ’wo’) 
grid on 

axis([0 24 -80 80]) 

title (’SEUS Zulu 18 Jam - Plot of TDE vs Time (Strip Chart’) 



y. LORAN-C Thesis 

y. 

y, seus20.m 

y. 
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% Data Set : Malone-SEUS 20 Jan97 

y. Program generates LORAN Time Difference Error (TDE) vs Time strip 
y, charts similar to the output of the current CALOC strip charts, 
y. These strip charts will be used to compare the output of the new 
y, control algorithm to determine the increase in control accuracy. 

y, Clear stored variables 
clear 

y. Load data set 
load seus20jn.dat; 

y, Remove bias from data. Original data has a range of 0-99 which 
y, equates to the actual range of -50 to +49. Each single step 
y, increase in the data equates to an actual 10 ns increase in the 
y. TDE, therefore the data is scaled by 10. 
seus20jn = (seus20jn - 50) ♦ 10; 

y. Load LORAN chain data vectors for Whiskey, Xray, Yainkee, aind Zulu. 

seus20jnW = seus20jn( ; , 1) ; 

seus20jnX = seus20jn( : ,2) ; 

seus20jnY = seus20jn( : ,8) ; 

seus20jnZ = seus20jn( : ,7) ; 

y. Determine the size of the data file. 

[num_rows, num_cols] = size(seus20jn) ; 

y. Determine number of data points after taking 7.5 minute averages. 
y« Since each data point represents 10 seconds in real time, there 
% are 45 data points in each 7.5 minute average. Using the Matlab 
y, floor command rounds down the result of the operation to get am 
y, integer value to use in reshaping the data for averaging. The 
y. second line generates the vector for plotting the data. 
data_cols = floor (num_rows/45) ; 
xaxis = [1 :num_rows/45] /8 ; 

y. Take 7.5 minute averages to generate data points for plotting 
y. similar to current strip charts . Averages are taken by reshaping 
y. the data into a matrix of dimensions data_cols x 45 and then tsdcing 
y. the mean in a column-wise fashion. This generates a vector of 
y. length data.cols of the 7.5 minute averages. 

avgW = meaii(reshape(seus20jnW(l : (45*data_cols)) , 45, data_cols) ) ; 
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avgX = mean(reshape(seus20jnX(l : (45*data_cols) ) , 45, data_cols) ) ; 
avgY = mean(reshape(seus20jnY(l : (45*data_cols) ) , 45, data_cols) ) ; 
avgZ = mean(reshape(seus20jnZ(l : (45*data_cols) ) , 45, data_cols)); 

*/. After plotting the raw data, the locations of the LPAs were 
*/. determined by comparing the generated plot with the actual strip 
*/. charts . The locations are placed below to generate LPA vectors . 
Ipaw = zeros (sized ;num_rows/45) ) ; 

Ipax = zeros (sized :num_rows/45) ) ; 

Ipay = zeros (sized :num_rows/45) ) ; 

Ipaz = zeros (sized :num_rows/45) ) ; 

Ipaw (30) = -20; 

Ipaw (44) = 20; Ipaw (89) = 20; Ipaw (104) = 20; 

Ipaw (123) = 20; Ipaw (160) = 20; 

Ipax (22) = -20; Ipax (68) = -20; Ipax (92) = -20; 

lpax(41) = 20; lpax(52) = 20; lpax(102) = 20; 

IpaxdOY) = 20; lpax(122) = 20; lpax(125) = 20; 

Ipax (142) = 20; Ipax (164) = 20; 

Ipay (7) = -20; Ipay (20) = -20; Ipay (35) = -20; 

Ipay (70) = -20; Ipay (87) = -20; 

lpay(96) = 20; Ipay(lOO) = 20; lpay(107) = 40; 

Ipaydll) = 20; lpay(122) = 20; lpay(124) = 40; 
lpay(142) = 20; lpay(164) = 20; 

Ipaz (28) = -20; Ipaz (49) = -20; 

Ipaz (87) = 20; Ipaz (98) = 20; Ipaz (104) = 20; 

Ipaz (110) = 20; Ipaz (122) = 20; Ipaz (126) = 20; 

Ipaz (128) = 20; Ipaz (131) = 20; Ipaz (140) = 20; 

Ipaz (164) = 20; 

y. Generate vector to plot LPA. Dividing by zero generates a 
y, vector of NaN so that only the LPA locations determined in 
y, the for loop below will be plotted, 
plotlpaw = zeros (size (Ipaw) )/0.0; 
plotlpax = zeros (size (Ipax) )/0.0; 



161 



plotlpay = zeros (size (lpay))/0.0; 
plotlpcLZ = zeros(size(lpaz))/0.0; 

% This for loop determines the locations of the LPAs by indexing 
% the capability of Matlab to index a vector based on its non-zero 
*/, values. The location of each LPA is transferred to the NaN 
7. vector generated above. Since Matlab does not plot NaN, the only 
% values that are plotted are the LPA locations, 
for index = l:data_cols 
if Ipaw (index) ~=0 

plotlpaw(index)=lpaw(index) ; 

end 

if Ipax (index) "=0 

plotlpax(index)=lpcLx(index) ; 

end 

if Ipay (index) ~=0 

plotlpay (index) =lpay (index) ; 

end 

if lpaz( index) ~=0 

plot lpaz( index) =lpaz (index) ; 

end 

end 

7, Generate the LORAN Strip Charts 
7. 

7. The plots are generated using subplot so that they appear with 
7. similar resolution to the existing strip charts. The plots 
7. include the TDE averages (avgY) , the cumulative TDE, and the 
7. locations of the LPAs input by CALOC. The initial value for 
7o the cumulative error is teiken from the actual strip charts from 
7o CALOC. The factor of 0.125 on the cumulative error is due to 
7 the frequency of the averages (i.e., 7.5 minutes is 1/8 of an 
7 hour) . The LPAs are plotted with a multiplicative factor of 3 
7 solely for ease of display, 
figure 

subplot (311) , plot(xaxis, avgW, ’w+’) 
xlabel(’Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, -40 + cumsum(avgW * 0.125), ’w-’) 
plot(xaxis, plot Ipaw * 3, ’wo’) 
grid on 

axis([0 24 -80 80]) 
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title (’SEUS Whiskey 20Jcin - Plot of TDE vs Time’) 



figure 

subplot (311) , plot(xaxis, avgX, ’w+’) 
xlabel(’Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, -60 + curasum(avgX * 0.125), ’w-’) 
plot(xaxis, plotlpax * 3, ’wo’) 
grid on 

axis([0 24 -80 80]) 

title (’SEUS Xray 20Jcin - Plot of TDE vs Time’) 
figure 

subplot (31 1) , plot(xaxis, avgY, ’w+’) 
xlabel(’Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, -50 + cumsum(avgY * 0.125), ’w-’) 
plot(xaxis, plotlpay * 3, ’wo’) 
grid on 

axis([0 24 -80 80]) 

title (’SEUS Yankee 20Jcin - Plot of TDE vs Time’) 
f igure 

subplot(311) , plot(xaxis, avgZ, ’w+’) 
xlabeK’Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, -92 + curasum(avgZ * 0.125), ’w-’) 
plot(xaxis, plotlpaz * 3, ’wo’) 
grid on 

axis([0 24 -80 80]) 

title (’SEUS Zulu 20Jan - Plot of TDE vs Time’) 



% LORAN-C Thesis 

7 . 

7, socuslS.m 

7. 

7o Data Set; Malone-SOCUS 18JcUi97 
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*/. Prograiin generates LORAN Time Difference Error (TDE) vs Time strip 
*/. charts similar to the output of the current CALOC strip charts. 

*/, These strip charts will be used to compare the output of the new 
*/, control algorithm to determine the increase in control accuracy. 

*/, Clear stored variables 
clear 

y. Load data set 
load socusl8j.dat; 

y, Remove bias from data. Original data has a range of 0-99 which 
y, equates to the actual range of -50 to +49. Each single step 
y, increase in the data equates to an actual 10 ns increase in the 
y, TDE, therefore the data is scaled by 10. 
socuslSj = (socuslSj - 50) * 10; 

y. Load LORAN chain data vectors for Victor, Whiskey, Xray, Yankee 
y, and Zulu. 

socuslSjV = socuslSj (:, 1) ; 
socuslSjW = socuslSj (: ,2) ; 
socuslSjX = socuslSj (: ,4) ; 
socuslSjY = socuslSj (: ,6) ; 
socuslSjZ = socuslSj (: ,9) ; 

y. Determine the size of the data file. 

[num.rows, nmn.cols] = size(socuslSj) ; 

y. Determine number of data points after taking 7.5 minute averages, 
y, Since each data point represents 10 seconds in real time, there 
y, are 45 data points in each 7.5 minute average. Using the Matlab 
y, floor command rounds down the result of the operation to get an 
y. integer value to use in reshaping the data for averaging. The 
y. second line generates the vector for plotting the data, 
data.cols = floor (num_rows/45) ; 
xaxis = [1 :num_rows/45]/S; 

y, Tadce 7.5 minute averages to generate data points for plotting 
y, similar to current strip charts. Averages are taken by reshaping 
y, the data into a matrix of dimensions data_cols x 45 and then taking 
y, the mean in a column-wise fashion. This generates a vector of 
y, length data_cols of the 7.5 minute averages. 
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avgV = mean (reshape (socuslSjV (1 : (45*data_cols) ) , 45, data_cols) ) ; 
avgW = mean(reshape(socusl8jW(l : (45*data_cols) ) , 45, data_cols)); 
avgX = mean(reshape(socusl8jX(l : (45*data_cols) ) , 45, data_cols)); 
avgY = mean(reshape(socusl8jY(l : (45*data_cols) ) , 45, data.cols)); 
avgZ = mean(reshape(socusl8jZ(l : (45*data_cols) ) , 45, data_cols)); 

% After plotting the raw data, the locations of the LPAs were 
7, determined by comparing the generated plot with the actual strip 
7 charts. The locations are placed below to generate LPA vectors. 
Ipav = zeros(size(l :num_rows/45) ) ; 

Ipaw = zeros (size (1 :num_rows/45) ) ; 

Ipax = zeros(size(l :num_rows/45)) ; 

Ipay = zeros (size (1 :num_rows/45) ) ; 

Ipaz = zeros(size(l :num_rows/45) ) ; 

lpav(94) = -20; 

Ipav (60) = 20; Ipav (132) = 20; Ipav (154) = 20; 

Ipaw (39) = -20; 

Ipaw (64) = 20; Ipaw (119) = 20; Ipaw (134) = 20; Ipaw (152) = 20; 

Ipax (24) = 20; 

Ipay (39) = -20; Ipay (132) = -20; 

Ipay (26) = 20; Ipay (30) = 20; 

Ipaz(lll) = -20; lpaz(113) = -20; lpaz(130) = -40; 

lpaz(12) = 20; lpaz(41) = 20; lpaz(65) = 20; lpaz(96) = 20; 

Ipaz (143) = 20; Ipaz (163) = 20; 

7 Generate vector to plot LPA. Dividing by zero generates a 

7 vector of NaN so that only the LPA locations determined in 

7 the for loop below will be plotted. 

plotlpav = zeros (size (Ipav) )/0.0; 

plotlpaw = zeros (size (Ipaw) )/0.0; 

plotlpax = zeros (size(lpax))/0.0; 

plotlpay = zeros (size (Ipay) )/0.0; 

plotlpaz = zeros(size(lpaz))/0.0; 
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y. This for loop determines the locations of the LPAs by indexing 
y. the capability of Matlab to index a vector based on its non-zero 
y, values . The location of each LPA is transferred to the NaN 
y, vector generated above. Since Matlab does not plot NaN, the only 
y, values that are plotted are the LPA locations, 
for index=l : data_cols 
if Ipav (index) ~=0 

plotlpav(index)=lpav(index) ; 

end 

if Ipaw ( index) ~=0 

plotlpaw(index)=lpaw(index) ; 

end 

if Ipax (index) ~=0 

plotlpax(index)=lpax(index) ; 

end 

if Ipay (index) ~=0 

plotlpay(index)=lpay(index) ; 

end 

if Ipaz (index) ~=0 

plotlpaz(index)=lpaz(index) ; 

end 

end 

y. Generate the LORAN Strip Charts 

y. 

y. The plots are generated using subplot so that they appear with 
y, similar resolution to the existing strip charts . The plots 
y, include the TDE averages (avgY) , the cumulative TDE, ctnd the 
y, locations of the LPAs input by CALOC. The initial value for 
y, the cumulative error is taken from the actual strip charts from 
y, CALOC. The factor of 0.125 on the cumulative error is due to 
y, the frequency of the averages (i.e., 7.5 minutes is 1/8 of an 
y, hour) . The LPAs are plotted with a multiplicative factor of 2 
y, or 3 solely for ease of display, 
figure 

subplot (311) , plot(xaxis, avgV, ’w+’) 
xlabel('Time (hours)'), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, -55 + cumsum(avgV * 0.125), 'w-') 
plot(xaxis, plotlpav * 3, 'wo') 
grid on 
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axis([0 24 -80 80]) 

title ('SOCUS Victor 18 Jan - Plot of IDE vs Time’) 



figure 

subplot (311) , plotCxaxis, avgW, ’w+’) 
xlabeK’Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, -17 + cumsum(avgW * 0.125), ’w-’) 
plot(xaxis, plotlpaw * 3, ’wo’) 
grid on 

axis([0 24 -80 80]) 

title (’SOCUS Whiskey 18 Jan - Plot of TDE vs Time’) 
figure 

subplot(311) , plot(xaxis, avgX, ’w+’) 
xlabeK’Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, 21 + cumsum(avgX * 0.125), ’w-’) 
plot(xaxis, plotlpax * 3, ’wo’) 
grid on 

axis([0 24 -80 80]) 

title (’SOCUS Xray 18 Jan - Plot of TDE vs Time’) 
figure 

subplot (311) , plot(xaxis, avgY, ’w+’) 
xlabeK’Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, 16 + cumsum(avgY * 0.125), ’w-’) 
plot(xaxis, plotlpay +3, ’wo’) 
grid on 

axis([0 24 -80 80]) 

title (’SOCUS Yankee 18 Jan - Plot of TDE vs Time’) 
figure 

subplot (311) , plot(xaxis, avgZ, ’w+’) 

XlabeK’Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, -20 + cumsum(avgZ * 0.125), ’w-’) 
plot(xaxis, plotlpaz *2, ’wo’) 
grid on 

axis([0 24 -80 80]) 
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title ('SOCUS Zulu 18 Jan - Plot of IDE vs Time') 



% LORAN-C Thesis 

% 

% socus20.m 

•/, 

% Data Set: Malone-SOCUS 20Jan97 

7, Program generates LORAN Time Difference Error (TDE) vs Time strip 
7, charts similar to the output of the current CALOC strip charts. 

7» These strip charts will be used to compare the output of the new 
7. control algorithm to determine the increase in control accuracy. 

7o Cleair stored variables 
clear 

7o Load data set 
load socus20j.dat; 

7o Remove bias from data. Original data has a range of 0-99 which 
7« equates to the actual range of -50 to +49. Each single step 
% increase in the data equates to an actual 10 ns increase in the 
7o TDE, therefore the data is scaled by 10. 
socus20j = (socus20j - 50) * 10; 

7. Load LORAN chain data vectors for Victor, Whiskey, Xray, Yankee 
7. and Zulu . 

socus20jV = socus20j ( : , 1) ; 
socus20jW = socus20j ( : ,2) ; 
socus20jX = socus20j ( : ,4) ; 
socus20jY = socus20j ( : ,6) ; 
socus20jZ = socus20j ( : ,9) ; 

7. Determine the size of the data file. 

[num_rows , num_cols] = size(socus20j) ; 

7. Determine number of data points after taking 7.5 minute averages. 
7. Since each data point represents 10 seconds in real time, there 
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% are 45 data points in each 7.5 minute average. Using the Matlab 
% floor command rounds down the result of the operation to get cin 
y, integer value to use in reshaping the data for averaging. The 
y, second line generates the vector for plotting the data, 
data.cols = floor (num_rows/45) ; 
xaxis = [l:num_rows/45]/8; 

*/o Take 7.5 minute averages to generate data points for plotting 
"/o similar to current strip charts. Averages are taken by reshaping 
y, the data into a matrix of dimensions data_cols x 45 aind then taking 
y. the mean in a column-wise fashion. This generates a vector of 
y, length data_cols of the 7.5 minute averages. 

avgV = mean(reshape(socus20jV(l: (45*data_cols)) , 45, data.cols)); 
avgW = mean(reshape(socus20jW(l; (45*data_cols)) , 45, data_cols)); 
avgX = mean(reshape(socus20jX(l; (45*data_cols)) , 45, data_cols)); 
avgY = mean(reshape (socus20jY(l : (45*data_cols) ) , 45, data_cols)); 
avgZ = mean(reshape(socus20jZ(l: (45*data_cols)) , 54, data_cols)); 

y, After plotting the raw data, the locations of the LPAs were 
y. determined by comparing the generated plot with the actual strip 
y, charts. The locations are placed below to generate LPA vectors. 
Ipav = zeros (size (1 : num_rows/45) ) ; 

Ipaw = zeros (size (1 : num_rows/45) ) ; 

Ipax = zeros (sized ;num_rows/45) ) ; 

Ipay = zeros (size (1 ;num_rows/45) ) ; 

IpcLZ = zeros (size (1 ;num_rows/45) ) ; 

Ipav (5) = 20; Ipav (133) = 20; 

lpaw(31) = -20; lpaw(69) = -20; 

lpaw(9) = 20; lpaw(14) = 20; lpaw(87) = 20; lpaw(143) = 20; 

Ipax (61) = 20; 
lpay(151) = -20; 

Ipay (64) = 20; Ipay (79) = 20; 

lpaz(105) = -20; lpaz(144) = -20; lpaz(158) = -20; 
lpaz(9) = 20; lpaz(77) = 20; lpaz(91) = 20; 
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*/. Generate vector to plot LPA. Dividing by zero generates a 

% vector of NaN so that only the LPA locations determined in 

y. the for loop below will be plotted. 

plotlpav = zeros (size (lpav))/0.0; 

plotlpaw = zeros (size (lpaw))/0.0; 

plotlpax = zeros (size (lpax))/0.0; 

plotlpay = zeros (size (lpay))/0.0; 

plotlpaz = zeros (size (lpaz))/0.0; 

y. This for loop determines the locations of the LPAs by indexing 
y. the capability of Matlab to index a vector based on its non-zero 
y. values . The location of each LPA is transferred to the NaN 
y, vector generated above. Since Matlab does not plot NaN, the only 
y, values that are plotted are the LPA locations, 
for index=l :data_cols 
if Ipav (index) "=0 

plotlpav(index)=lpav(index) ; 

end 

if Ipaw (index) ~=0 

plotlpaw(index)=lpaw(index) ; 

end 

if Ipax (index) "=0 

plotlpeix (index) =lpax (index) ; 

end 

if Ipay ( index) ~=0 

plotlpay(index)=lpay(index) ; 

end 

if Ipciz (index) ~=0 

plotlpaz(index)=lpaz(index) ; 

end 

end 

y. Generate the LORAN Strip Charts 

y. 

y. The plots cire generated using subplot so that they appear with 
y, similar resolution to the existing strip charts. The plots 
y. include the TDE averages (avgY) , the cumulative TDE, and the 
y, locations of the LPAs input by CALOC. The initial value for 
y, the cumulative error is taken from the actual strip charts from 
y, CALOC. The factor of 0.125 on the cumulative error is due to 
y, the frequency of the averages (i.e., 7.5 minutes is 1/8 of an 
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y, hour) . The LPAs are plotted with a multiplicative factor of 2 

y. or 3 solely for ease of display. 

figure 

subplotOll) , plot(xaxis, avgV, ’w+’) 
xlabelC’Time (hours)’), ylabelC’TDE (ns)’) 
hold on 

plot(xaxis, -42 + cumsum(avgV * 0.125), ’w-’) 
plot(xaxis, plotlpav * 3, ’wo’) 
grid on 

axis([0 24 -80 80]) 

title(’S0CUS Victor 20 Jan - Plot of TDE vs Time’) 
figure 

subplot (311) , plot(xaxis, avgW, ’w+’) 
xlabeK’Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, -20 + cumsum(avgW * 0.125), ’w-’) 
plot(xaxis, plotlpaw * 3, ’wo’) 
grid on 

axis([0 24 -80 80]) 

title (’SOCUS Whiskey 20Jan - Plot of TDE vs Time’) 
figure 

subplot (311) , plot(xaxis, avgX, ’w+’) 
xlabel(’Time (hours)’), ylabel(’TDE (ns)’) 
hold on 

plot(xaxis, 20 + ciunsiun(avgX * 0.125), ’w-’) 
plot(xaxis, plotlpax * 3, ’wo’) 
grid on 

axis([0 24 -80 80]) 

title (’SOCUS Xray 20 Jan - Plot of TDE vs Time’) 
figure 

subplot (31 1) , plot(xaxis, avgY, ’w+’) 
xlabeK’Time (hours)’), ylabel(’TDE (ns)’) 
hold on 

plot(xaxis, -42 + cumsum(avgY * 0.125), ’w-’) 
plot(xaxis, plotlpay * 3, ’wo’) 
grid on 

axis([0 24 -80 80]) 

title(’ SOCUS Yankee 20Jan - Plot of TDE vs Time’) 
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figure 

subplot (311) , plot(xaxis, avgZ, ’w+') 
xlabeK’Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, 4 + cumsum(avgZ * 0.125), ’w-’) 
plot(xaxis, plotlpaz * 2, ’wo’) 
grid on 

axis ( [0 24 -80 80] ) 

title (’SOCUS Zulu 20 Jan - Plot of TDE vs Time’) 



•/. LORAN-C Thesis 

•/. 

'/, nocus4.m 

7 . 

7. Data Set : Middletown-NOCUS 4May97 

7o Program generates LORAN Time Difference Error (TDE) vs Time strip 
7. charts similar to the output of the current CALOC strip charts . 

7o These strip charts will be used to compare the output of the new 
7. control algorithm to determine the increase in control accuracy. 

7. Clear stored variables 
clear 

7. Load data set 
load nocus4my.dat; 

7. Remove bias from data. Original data has a range of 0-99 which 
7« equates to the actual range of -50 to +49. Each single step 
7. increase in the data equates to am actual 10 ns increase in the 
7 TDE, therefore the data is scaled by 10. 
nocus4my = (nocus4my - 50) * 10; 

7o Load LORAN chain data vectors for Whiskey, Xray and Yankee. 
nocus4myW = nocus4my( : , 10) ; 
nocus4myX = nocus4my( : , 12) ; 
nocus4myY = nocus4my( : ,3) ; 
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7. Determine the size of the data file. 

[num_rows, num_cols] = size (nocus4my) ; 

7. Determine number of data points after taking 7.5 minute averages. 

7. Since each data point represents 10 seconds in real time, there 
7 are 45 data points in each 7.5 minute average. Using the Matlab 
7o floor command rounds down the result of the operation to get an 
7. integer value to use in reshaping the data for averaging. The 
7. second line generates the vector for plotting the data, 
data.cols = floor (num_rows/45) ; 
xaxis = [1 :nvim_rows/45]/8; 

7o Take 7.5 minute averages to generate data points for plotting 
7o similar to current strip charts. Averages are taken by reshaping 
7o the data into a matrix of dimensions data_cols x 45 and then talcing 
7o the mean in a column-wise fashion. This generates a vector of 
7. length data_cols of the 7.5 minute averages. 

avgW = meaui(reshape(nocus4myW(l : (45*data_cols) ) , 45, data_cols)); 
avgX = mecin(reshape(nocus4myX(l: (45*data_cols)) , 45, data_cols)); 
avgY = meein(reshape(nocus4myY(l : (45*data_cols) ) , 45, data_cols)); 

7. After plotting the raw data, the locations of the LPAs were 
7o determined by comparing the generated plot with the actual strip 
7o charts. The locations are placed below to generate LPA vectors. 
Ipaw = zeros(size(l:num_rows/45)) : 

Ipax = zeros (size (1 :num_rows/45) ) ; 

Ipay = zeros(size(l :num_rows/45) ) ; 

Ipaw (108) = -20; Ipaw (133) = -20; 
lpaw(157) = -20; lpaw(178) = -20; 

Ipaw (13) = 40; Ipaw (24) = 20; Ipaw (79) = 20; 

Ipax (145) = -20; Ipax (157) = -20; 

Ipax (48) = 20; Ipax (178) = 20; 

Ipay (126) = -20; lpay(146) = -20; 

Ipay (158) = -20; Ipay (178) = -20; 

Ipay (39) = 20; Ipay (55) = 20; Ipay (60) = 40; 
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*/, Generate vector to plot LPA. Dividing by zero generates a 

7. vector of NaN so that only the LPA locations determined in 

7, the for loop below will be plotted. 

plotlpaw = zeros (size (lpaw))/0.0; 

plotlpcLX = zeros (size (lpcLx))/0.0: 

plotlpay = zeros (size (lpay))/0.0; 

7» This for loop determines the locations of the LPAs by indexing 
7o the capability of Matlab to index a vector based on its non-zero 
7. values. The location of each LPA is transferred to the NaN 
7. vector generated above. Since Matlab does not plot NaN, the only 
7« values that are plotted are the LPA locations, 
for index = l;data_cols 
if lpaw( index) ~=0 

plotlpaw(index)=lpaw(index) ; 

end 

if Ipax ( index) "=0 

plotlpax( index) =lpax( index) ; 

end 

if Ipay (index) "=0 

plotlpay ( index) =lpay( index) ; 

end 

end 

7 Generate the LORAN Strip Charts 

7. 

7. The plots are generated using subplot so that they appear with 
7. similar resolution to the existing strip charts. The plots 
7. include the TDE averages (avgY), the cumulative TDE, and the 
7. locations of the LPAs input by CALOC. The initial value for 
7. the cumulative error is taJcen from the actual strip charts from 
7. CALOC. The factor of 0.125 on the cumulative error is due to 
7. the frequency of the averages (i.e., 7.5 minutes is 1/8 of an 
7. hour) . The LPAs are plotted with a multiplicative factor of 3 
7. solely for ease of display, 
figure 

subplot (311) , plot(xaxis, avgW, 'w+’) 
xlabeK'Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, -35 + cumsum(avgW =1= 0.125), ’w-’) 
plot(xaxis, plotlpaw * 3, ’wo’) 
grid on 
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axis([0 24 -80 80]) 

title (’NOCUS Whiskey 4May - Plot of TDE vs Time’) 
figure 

subplotOll) , plotCxaxis, avgX, ’w+’) 
xlabeK’Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, 30 + cumsum(avgX * 0.125), ’w-’) 
plot(xaxis, plotlpax * 3, ’wo’) 
grid on 

axis([0 24 -80 80]) 

title (’NOCUS Xray 4May - Plot of TDE vs Time’) 
figure 

subplot(311) , plot(xaxis, avgY, ’w+’) 
XlabeK’Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, 25 + cumsum(avgY * 0.125), ’w-’) 
plot(xaxis, plotlpay * 3, ’wo’) 
grid on 

axis([0 24 -80 80]) 

title (’NOCUS Yankee 4May - Plot of TDE vs Time’) 



y. LORAN-C Thesis 

y. 

y. nocus5.m 

y. 

y. Data Set; Middletown-NOCUS 5May97 

y. Program generates LORAN Time Difference Error (TDE) vs Time strip 
% charts similar to the output of the current CALOC strip charts, 
y, These strip charts will be used to compare the output of the new 
y control algorithm to determine the increase in control accuracy. 

y Clear stored variables 
clear 

y Load data set 
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load nocus5my.dat; 



*/. Remove bias from data. Original data has a rainge of 0-99 which 
*/. equates to the actual range of -50 to +49. Each single step 
*/. increase in the data equates to an actual 10 ns increase in the 
*/, IDE, therefore the data is scaled by 10. 
nocus5my = (nocusSmy - 50) * 10; 

% Load LORAN chain data vectors for Whiskey, Xray and Yankee. 
nocus5myW = nocus5my( : , 10) ; 
nocus5myX = nocus5my( : , 12) ; 
nocus5myY = nocus5my( : ,3) ; 

’/, Determine the size of the data file. 

[num_rows, num_cols] = size(nocus5my) ; 

y, Determine number of data points after taking 7.5 minute averages, 
y, Since each data point represents 10 seconds in real time, there 
y, are 45 data points in each 7.5 minute average. Using the Matlab 
y, floor command rounds down the result of the operation to get an 
y, integer value to use in reshaping the data for averaging. The 
y, second line generates the vector for plotting the data. 
data_cols = floor( num_rows/45 ) ; 
xaxis = [1 :num_rows/45]/8; 

y. Take 7.5 minute averages to generate data points for plotting 
y, similar to current strip charts. Averages are taken by reshaping 
y, the data into a matrix of dimensions data_cols x 45 and then taking 
y, the mean in a column-wise fashion. This generates a vector of 
y, length data_cols of the 7 . 5 minute averages . 

avgW = mean(reshape(nocus5myW(l : (45*data_cols) ) , 45, data_cols)) ; 
avgX = mecLn(reshape(nocus5myX(l : (45*data_cols) ) , 45, data_cols)); 
avgY = mean(reshape(nocus5myY(l : (45*data_cols) ) , 45, data_cols)); 

y. After plotting the raw data, the locations of the LPAs were 
y, determined by comparing the generated plot with the actual strip 
y. charts. The locations are placed below to generate LPA vectors. 
Ipaw = zerosCsized :num_rows/45) ) ; 

Ipax = zeros (sized : num_rows/45) ) ; 

Ipay = zerosCsized :num_rows/45) ) ; 

lpaw(147) = -20; 
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lpax(105) = -20; 



lpax(12) = 20; 

lpay(115) = -20; lpay(147) = -40; 

lpay(13) = 20; lpay(30) = 20; 
lpay(50) = 20; lpay(74) = 20; 

*/, Generate vector to plot LPA. Dividing by zero generates a 

7, vector of NaN so that only the LPA locations determined in 

7, the for loop below will be plotted. 

plotlpaw = zeros (size (Ipaw) )/0.0; 

plotlpax = zeros(sizedpax) )/0.0; 

plotlpay = zeros (size (Ipay) )/0.0; 

% This for loop determines the locations of the LPAs by indexing 
7. the capability of Matlab to index a vector based on its non-zero 
7 values. The location of each LPA is transferred to the NaN 
7. vector generated above. Since Matlab does not plot NaN, the only 
7. values that are plotted are the LPA locations . 
for index = l:data_cols 
if Ipaw (index) "=0 

plotlpaw(index)=lpaw(index) ; 

end 

if Ipax (index) ~=0 

plotlpax(index)=lpax(index) ; 

end 

if Ipay (index) ~=0 

plotlpay ( index) =lpay (index) ; 

end 

end 

7« Generate the LORAN Strip Charts 
7. 

7o The plots are generated using subplot so that they appear with 
7. similar resolution to the existing strip charts. The plots 
7» include the TDE averages (avgY) , the cumulative TDE, and the 
7. locations of the LPAs input by CALOC. The initial value for 
7. the cumulative error is taken from the actual strip charts from 
7. CALOC. The factor of 0.125 on the cumulative error is due to 
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7. the frequency of the averages (i.e., 7.5 minutes is 1/8 of an 
7o hour) . The LPAs are plotted with a multiplicative factor of 3 
7. solely for ease of display, 
figure 

subplot(311) , plotCxaxis, avgW, ’w+’) 
xlabeK’Time (hours)'), ylabeK'TDE (ns)') 
hold on 

plot(xaxis, -15 + curasum(avgW * 0.125), 'w-') 
plot(xaxis, plotlpaw * 3, 'wo') 
grid on 

axis([0 24 -80 80]) 

title ('NOCUS Whiskey 5May - Plot of TDE vs Time') 
figure 

subplot (311) , plot(xaxis, avgX, 'w+') 
xlabeK'Time (hours)'), ylabeK'TDE (ns)') 
hold on 

plot(xaxis, -25 + cums\im(avgX * 0.125), 'w-') 
plot(xaxis, plotlpax * 3, 'wo') 
grid on 

axis([0 24 -80 80]) 

title ('NOCUS Xray 5May - Plot of TDE vs Time') 
figure 

subplot (311) , plot(xaxis, avgY, 'w+') 
xlabeK'Time (hours)'), ylabeK'TDE (ns)') 
hold on 

plot(xaxis, 5 + cumsum(avgY * 0.125), 'w-') 
plot(xaxis, plotlpay * 3, 'wo') 
grid on 

axis([0 24 -80 80]) 

title ('NOCUS Ycinkee 5May - Plot of TDE vs Time') 



7. LORAN-C Thesis 
7. 

7. uswc4.m 
7. 
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7. Data Set: USWC 4May97 

7o Program generates LORAN Time Difference Error (TDE) vs Time strip 
7, charts similar to the output of the current CALOC strip charts. 

7o These strip charts will be used to compare the output of the new 
7o control algorithm to determine the increase in control accuracy. 

7o Cleeir stored variables 
clear 

7o Load data set 
load uswc4my . dat ; 

7o Remove bias from data. Original data has a range of 0-99 which 
7o equates to the actual range of -50 to +49. Each single step 
7. increase in the data equates to ein actual 10 ns increase in the 
7. TDE, therefore the data is scaled by 10. 
uswc4my = (uswc4my - 50) * 10; 

7o Load LORAN chain data vectors for Whiskey, Xray and Yankee. 
uswc4myW = uswc4my( : , 1) ; 
uswc4myX = uswc4my( : , 10) ; 
uswc4myY = uswc4my( : , 12) ; 

7. Determine the size of the data file. 

[num_rows, num_cols] = size (uswc4my) ; 

7. Determine number of data points after taking 7.5 minute averages. 

7. Since each data point represents 10 seconds in real time, there 
7. are 45 data points in each 7.5 minute average. Using the Matlab 
7. floor command rounds down the result of the operation to get cin 
7. integer value to use in reshaping the data for averaging. The 
7. second line generates the vector for plotting the data. 
data_cols = floor ( num_rows/45 ); 
xaxis = [1 :num_rows/45]/8; 

7. Take 7.5 minute averages to generate data points for plotting 
7. similar to current strip charts. Averages are taken by reshaping 
7o the data into a matrix of dimensions data_cols x 45 and then talking 
7o the mean in a column-wise fashion. This generates a vector of 
7o length data_cols of the 7.5 minute averages. 

avgW = meau(reshape(uswc4myW(l: (45*data_cols)) , 45, data.cols)); 
avgX = mecin(reshape(uswc4myX(l: (45*data_cols)) , 45, data_cols)); 
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avgY = meaji(reshape (uswc4myY(l ; (45*data_cols) ) , 45, data_cols)); 



y. After plotting the raw data, the locations of the LPAs were 
y, determined by comparing the generated plot with the actual strip 
y, charts. The locations are placed below to generate LPA vectors. 
Ipaw = zeros (size (1 :num_rows/45)) 

Ipax = zeros (size (1 :num_rows/45)) 

Ipay = zeros (size (l:num_rows/45)) 

Ipaw (5) = -20; Ipaw (45) = -20; 

Ipaw (50) = -20; Ipaw (84) = -20; 

Ipaw (94) = 20; Ipaw (104) = 20; 

lpauc(69) = -20; 

lpax(34) = 20; lpax(50) = 20; lpax(154) = 20; 

Ipay (84) = -20; 

Ipay (40) = 20; Ipay (55) = 20; 

y. Generate vector to plot LPA. Dividing by zero generates a 

y, vector of NaN so that only the LPA locations determined in 

y. the for loop below will be plotted. 

plotlpaw = zeros (size (Ipaw) )/0.0; 

plotlpax = zeros (size (Ipax) )/0.0; 

plotlpay = zeros (size (Ipay) )/0.0; 

y. This for loop determines the locations of the LPAs by indexing 
y, the capability of Matlab to index a vector based on its non-zero 
y, values. The location of each LPA is trainsf erred to the NaN 
y, vector generated above. Since Matlab does not plot NaN, the only 
y, values that are plotted are the LPA locations, 
for index = l:data_cols 
if Ipaw ( index) ~=0 

plotlpaw(index)=lpaw(index) ; 

end 

if Ipax(index) ~=0 

plotlpax ( index) =lpax (index) ; 

end 

if Ipay (index) ~=0 



180 



plot lpay( index) =lpay( index) ; 



end 

end 

*/. Generate the LORAN Strip Charts 

7. 

7, The plots are generated using subplot so that they appear with 
7. similar resolution to the existing strip charts . The plots 
7. include the TDE averages (avgY) , the cumulative TDE, cind the 
7. locations of the LPAs input by CALOC. The initial value for 
7. the cumulative error is taken from the actual strip charts from 
7. CALOC. The factor of 0.125 on the cumulative error is due to 
7. the frequency of the averages (i.e., 7.5 minutes is 1/8 of cin 
7o hour) . The LPAs are plotted with a multiplicative factor of 3 
7. solely for ease of display, 
figure 

subplot (311) , plotCxaxis, avgW, ’w+’) 
xlabeK’Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, 80 + cumsum(avgW * 0.125), ’w-’) 
plot(xaxis, plotlpaw * 3, ’wo’) 
grid on 

axis([0 24 -80 80]) 

title(’USWC Whiskey 4May - Plot of TDE vs Time’) 
figure 

subplot(311) , plot(xaxis, avgX, ’w+’) 
xlabel(’Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, -30 + cumsum(avgX * 0.125), ’w-’) 
plot(xaxis, plotlpax * 3, ’wo’) 
grid on 

axis([0 24 -80 80]) 

title(’USWC Xray 4May - Plot of TDE vs Time’) 
figure 

subplot (311) , plot(xaxis, avgY, ’w+’) 
xlabeK’Time (hours)’), ylabel(’TDE (ns)’) 
hold on 

plot(xaxis, -40 + cumsum(avgY * 0.125), ’w-’) 
plot(xaxis, plotlpay * 3, ’wo’) 
grid on 
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axis([0 24 -80 80]) 

title (’USWC Yankee 4May - Plot of IDE vs Time') 



•/. LORAN-C Thesis 

y. 

y, uswcS.m 

y. 

y. Data Set : Middletown-USWC 5May97 

y. Program generates LORAN Time Difference Error (TDE) vs Time strip 
y, charts similair to the output of the current CALOC strip charts, 
y. These strip charts will be used to compare the output of the new 
y. control algorithm to determine the increase in control accuracy. 

y. Clear stored variables 
clear 

y. Load data set 
load uswc5my.dat; 

y. Remove bias from data. Original data has a range of 0-99 which 
y, equates to the actual range of -50 to +49. Each single step 
y. increase in the data equates to an actual 10 ns increase in the 
y, TDE, therefore the data is scaled by 10. 
uswcSmy = (uswcSmy - 50) * 10; 

y. Load LORAN chain data vectors for Whiskey, Xray and Yankee. 
uswc5myW = uswc5my ( : , 1) ; 
uswc5rayX = uswc5my ( : , 10) ; 
uswc5myY = uswc5my ( : , 12) ; 

y. Determine the size of the data file. 

[n\im_rows, num_cols] = size(uswc5my) ; 

y. Determine number of data points after taking 7.5 minute averages, 
y. Since each data point represents 10 seconds in real time, there 
y are 45 data points in each 7.5 minute average. Using the Matlab 
y. floor command rounds down the result of the operation to get an 



182 



y, integer value to use in reshaping the data for averaging. 
7. second line generates the vector for plotting the data, 
data.cols = floor ( nuin_rows/45 ); 
xaxis = [1 ;num_rows/45] /8; 



The 



y. Take 7.5 minute averages to generate data points for plotting 
y, similar to current strip charts. Averages are taken by reshaping 
y, the data into a matrix of dimensions data_cols x 45 and then taking 
y. the meain in a column-wise fashion. This generates a vector of 
y. length data_cols of the 7.5 minute averages. 
avgW = meaii(reshape(uswc5myW(l ; (45*data_cols)) , 45, data_cols)) 
avgX = mecLn(reshape (uswc5myX(l : (45*data_cols) ) , 45, data_cols)) 
avgY = mecin (reshape (uswc5myY( 1 ; (45*data_cols) ) , 45, data_cols)) 

7, After plotting the raw data, the locations of the LPAs were 
y. determined by comparing the generated plot with the actual strip 
y. charts. The locations are placed below to generate LPA vectors. 
Ipaw = zeros (size (1 ;num_rows/45) ) ; 

Ipax = zeros(size(l :num_rows/45) ) ; 

Ipay = zeros (size (1 :num_rows/45) ) ; 

lpaw(45) = -20; lpaw(61) = -20; 

Ipaw (71) = -20; 

Ipaw (94) = 20; Ipaw (121) = 20; 

Ipax (75) = -40; 

Ipax (33) = 20; lpax(81) = 20; lpax(94) = 20; 

Ipay (75) = -40; 

Ipay (48) = 20; lpay(94) = 20; Ipay (157) = 20; 

y Generate vector to plot LPA. Dividing by zero generates a 
% vector of NaN so that only the LPA locations determined in 
y the for loop below will be plotted, 
plotlpaw = zeros (size (Ipaw) )/0.0; 
plotlpax = zeros (size (Ipax) )/0.0; 
plotlpay = zeros (size (Ipay) )/0.0; 

y This for loop determines the locations of the LPAs by indexing 
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y, the capability of Matlab to index a vector based on its non-zero 
y, values. The location of each LPA is treinsf erred to the NaN 
y, vector generated above. Since Matlab does not plot NaN, the only 
y, values that are plotted are the LPA locations, 
for index = l:data_cols 
if Ipaw(index) ~=0 

plotlpaw(index)=lpaw(index) ; 

end 

if Ipax (index) ~=0 

plotlpax(index)=lpax(index) ; 

end 

if Ipay (index) ~=0 

plotlpay (index) =lpay (index) ; 

end 

end 

y. Generate the LORAN Strip Charts 

y. 

y. The plots are generated using subplot so that they appear with 
y. similar resolution to the existing strip charts. The plots 
y, include the TDK averages (avgY) , the cumulative TDE, aind the 
% locations of the LPAs input by CALOC. The initial value for 
y, the cumulative error is taken from the actual strip charts from 
y, CALOC. The factor of 0.125 on the cumulative error is due to 
y, the frequency of the averages (i.e., 7.5 minutes is 1/8 of an 
y hour) . The LPAs are plotted with a multiplicative factor of 2 
y. or 3 solely for ease of display, 
figure 

subplot (311) , plot(xaxis, avgW, ’w+’) 
xlabeK’Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, 60 + cumsum(avgW * 0.125), ’w-’) 
plot(xaxis, plotlpaw +3, ’wo’) 
grid on 

axis([0 24 -80 80]) 

title (’USWC Whiskey 5May - Plot of TDE vs Time’) 
figure 

subplot (311) , plot(xaxis, avgX, ’w+’) 
xlabeK’Time (hours)’), ylabel(’TDE (ns)’) 
hold on 

plot(xaxis, 25 + cumsum(avgX * 0.125), ’w-’) 
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plotCxaxis, plotlpax * 2 , 'wo’) 
grid on 

axis([0 24 -80 80]) 

title (’USWC Xray 5May - Plot of TDE vs Time’) 
figure 

subplotOll) , plot(xaxis, avgY, ’w+’) 
xlabeK’Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, 25 + cumsum(avgY * 0.125), ’w-’) 
plot(xaxis, plotlpay * 2 , ’wo’) 
grid on 

axis([0 24 -80 80]) 

title (’USWC Yankee 5May - Plot of TDE vs Time’) 



•/. LORAN-C Thesis 

•/. 

*/, kodiaJczulu.m 

•/. 

’/. Data Set : KodiaJc Zulu 

7. Program generates LORAN Time Difference Error (TDE) vs Time strip 
7, charts similar to the output of the current CALOC strip charts. 

7o These strip charts will be used to compare the output of the new 
7« control algorithm to determine the increase in control accuracy. 

7. Clear stored variables 
clear 

7. Load data set 
load kodiaJc.dat; 

7. Remove bias from data. Original data has a range of 0-99 which 
7. equates to the actual range of -50 to +49. Each single step 
7. increase in the data equates to an actual 10 ns increase in the 
7. TDE, therefore the data is scaled by 10. 
kodiak = (kodieik - 50) * 10; 
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y. Load individual LORAN chain data vector for Zulu. 
kodieikZ = kodiak(:,l); 

y. Determine the size of the data file. 

[num_rows, num_cols] = size(kodicik) ; 

y, Determine number of data points after tciking 7.5 minute averages, 
y. Since each data point represents 10 seconds in real time, there 
y, are 45 data points in each 7.5 minute average. Using the Matlab 
y, floor command rounds down the result of the operation to get an 
% integer value to use in reshaping the data for averaging. The 
y. second line generates the vector for plotting the data. 
data_cols = floor( num_rows/45 ); 
xaxis = [1 :num_rows/45] /8; 

y. Take 7.5 minute averages to generate data points for plotting 
% similar to current strip charts . Averages are taiken by reshaping 
y, the data into a matrix of dimensions data_cols x 45 and then taking 
y, the mean in a column-wise fashion. This generates a vector of 
*/o length data_cols of the 7.5 minute averages. 
avgZ = mecui(reshape (kodicikZd : (45*data_cols) ) , 45, data_cols)); 

y. After plotting the raw data, the locations of the LPAs were 
y« determined by comparing the generated plot with the actual strip 
y, charts . The locations are placed below to generate LPA vectors . 
Ipaz = zeros (sized :data_cols) ) ; 

lpaz(37) = -20; Ipaz (82) = -20; 

Ipaz (9) = 20; Ipaz (46) = 20; 

y. Generate vector to plot LPA. Dividing by zero generates a 
y. vector of NaN so that only the LPA locations determined in 
y. the for loop below will be plotted, 
plotlpaz = zeros (size (Ipaz) )/0.0; 

y. This for loop determines the locations of the LPAs by indexing 
y, the capability of Matlab to index a vector based on its non-zero 
y, values. The location of each LPA is transferred to the NaN 
y, vector generated above. Since Matlab does not plot NaN, the only 
y, values that are plotted are the LPA locations, 
for index=l : data_cols 
if Ipaz (index) ~=0 
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plotlpaz (index) =lpaz (index) ; 



end 

end 

% Generate the LORAN Strip Charts 

7 . 

7. The plots are generated using subplot so that they appear with 
7. similar resolution to the existing strip charts . The plots 
7. include the TDE averages (avgY) , the cumulative TDE, and the 
7. locations of the LPAs input by CALOC. The initial value for 
7o the cumulative error is taken from the actual strip charts from 
7. CALOC. The factor of 0.125 on the cumulative error is due to 
7. the frequency of the averages (i.e., 7.5 minutes is 1/8 of an 
7o hour) . The LPAs are plotted with a multiplicative factor of 3 
7. solely for ease of display, 
subplot (311) , plot(xaxis, avgZ, ’w+’) 
xlabeK'Time (hours)’), ylabeK’TDE (ns)’) 
hold on 

plot(xaxis, -15 + cumsum(avgZ * 0.125), ’w-’) 

plot(xaxis, plotlpaz * 3, ’wo’) 

grid on, axis([0 24 -80 80]) 

title( ’Kodicik Zulu - Plot of TDE vs Time’) 
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7. PARSE TITLE FUNCTION 



function graph_title = parse_title (string) 

% Loran-C Thesis 

1 

7o PARSE_TITLE This function parses the data file name string cind 
*/. generates a string for the title of the output graph. 

*/, Get Lorcin-C chain name 
if strcmp(string(l :5) , ’seus_') 
chain = ’SEUS’; 

elseif strcmp (string (1 : 5) , ’neus_’) 
chain = ’NEUS’; 

elseif strcmp(string(l : 5) , ’socus’) 
chain = ’SOCUS’; 

elseif strcmp (string( 1 : 5) , ’nocus’) 
chain = ’NOCUS’; 

elseif strcrap(string(l : 5) , ’uswc_’) 
chain = ’USWC'; 

elseif strcmp(string(l :5) , ’kodik’) 
chain = ’KODIAK’; 

end 

y. Get date 

date = string(6 ;7) ; 

y. Get month 

if strcmp (string(8 : 9) , ’jn’) 
month = ’January’ ; 
elseif strcmp(string(8:9) , ’fb’) 
month = ’February’; 
elseif strcmp (string (8: 9) , ’mr’) 
month = ’March’; 

elseif strcmp(string(8:9) , ’ar’) 
month = ’April’; 

elseif strcmp (string (8: 9) , ’my’) 
month = ’May’; 

elseif strcmp (string (8; 9) , ’ju’) 
month = ’ June ’ ; 
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elseif strcmp(string(8:9) , ’jlO 
month = ’July’; 

elseif strcmp(string(8 : 9) , ’ag’) 
month = ’ August ’ ; 
elseif strcmp(string(8 : 9) , ’sp’) 
month = ’September’; 
elseif strcmp(string(8 ; 9) , ’ot’) 
month = ’October’; 
elseif strcmp(string(8 : 9) , ’nv’) 
month = ’November’; 
elseif strcmp(string(8:9) , ’dc’) 
month = ’December’; 

end 

% Get Loran-C station designator 
if strcmp(stringdO) , ’V’) 
station = ’Victor’; 
elseif strcmp(stringdO) , ’W’) 
station = ’Whiskey’; 
elseif strcmp(stringdO) , ’X’) 
station = ’Xray’; 
elseif strcmp(string(10) , ’Y’) 
station = ’Yankee’; 
elseif strcmp(stringdO) , ’Z’) 
station = ’Zulu’; 

end 

graph_title = [chain, ’ ’, station, ’ (’, date, ’ 



’ , month , ’ ) ’ ] ; 
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